import os
import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

import warnings
from astropy.coordinates import Distance
import astropy.constants as const
from astropy.time import Time
from astropy.io import fits


import warnings
import os
import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

import warnings
from astropy.coordinates import Distance
import astropy.constants as const
from astropy.time import Time
from astropy.io import fits
class curve_fit_mcmc():
    def __init__(self,f,x_data,y_data,y_err=[],bounds=[],
                 p0=[],
                 N_pool=4,seed=42,nwalkers=32,
                 nsamplers=5000,nburn=1500,save_chain=True,
                 prior=None,
                 ylimit=np.array([]),args=None
                 ):
        """
        f: 拟合的函数
        x_data: 数据x
        y_data: 数据y
        y_err: 数据y的误差
        bounds: 参数的边界，默认为空，即不限制参数的取值范围[-inf,inf]
        p0: 初始参数
        N_pool: 并行计算的进程数
        seed: 随机种子
        nwalker: 每个walker的参数个数
        nsamplers: 采样步数
        nburn: 预热步数
        save_chain: 是否保存采样结果
        ylimit :numpy array. 记录数据是上限(>0)/下限(<0)/实际探测值(0).
        prior: 先验方程，规定了函数之间内部的一些关系。

        """
        import inspect
        self.func=f
        self.x_data=x_data
        self.y_data=y_data
        self.bounds=bounds
        self.y_err =y_err
        self.args=args
        self.N_pool=N_pool
        self.p0 = np.array(p0)
        self.seed=seed
        self.ndim= len(p0)
        self.nwalker=nwalkers
        self.nsamplers=nsamplers
        self.nburn=nburn
        self.save_chain=save_chain
        self.ylimit=ylimit if np.any(ylimit) else np.zeros(len(y_data))
        self.func_name=self.func.__name__
        self.prior=self.prior_func if prior==None else prior
        self.check_prior=False if prior==None else True
        self.check_limit=True if np.any(np.where(self.ylimit!=0)[0]) else False


        self.upperlimit_arg=np.where(self.ylimit>0)[0]
        self.lowerlimit_arg=np.where(self.ylimit<0)[0]
        self.real_data_arg =np.where(self.ylimit==0)[0]
        self.x_data_real=self.x_data[self.real_data_arg]
        self.y_data_real=self.y_data[self.real_data_arg]
        self.y_err_real =self.y_err [self.real_data_arg]
        self.check_y_err()
        # self.run_mcmc()
    def prior_func(self,params): 
        return 0
        
    def check_limit_data(self,y_model):
        if np.any(self.upperlimit_arg):
            if np.any(y_model[self.upperlimit_arg]>self.y_data[self.upperlimit_arg]):
                return -np.inf
        if np.any(self.lowerlimit_arg):

            if np.any(y_model[self.lowerlimit_arg]<self.y_data[self.lowerlimit_arg]):

                return -np.inf
        return 0
    
    def check_bound(self):
        if not any(list(self.bounds)):
            self.bounds=[tuple(self.ndim*[-np.inf]),
                         tuple(self.ndim*[ np.inf])]
    def check_y_err(self):
        self.y_err=np.abs(self.y_err)
        if np.any(self.y_err == 0):
            self.y_err[self.y_err==0]=1e-5
            warnings.warn("y_err has 0, and which have set to 1e-5.")
        if np.any(np.isnan(self.y_err)):
            self.y_err[np.isnan(self.y_err)]=1
            warnings.warn("y_err has None, and which have set to 1.")
        if np.any(np.isposinf(self.y_err)):
            self.y_err[np.isposinf(self.y_err)]=10
            warnings.warn("y_err has inf, and which have set to 10.")

    def check_bound(self,params):
        for i,p in enumerate(params):
            if not self.bounds[0][i]<p<self.bounds[1][i]:
                return -np.inf
        return 0
    def log_prob(self,params,*args):
        if np.isneginf(self.check_bound(params)): return -np.inf
        if self.check_prior: 
            
            if np.isneginf(self.prior_func (params)):
                return -np.inf
        y_model=self.func(params,self.x_data,args)
        if self.check_limit:
            
            if np.isneginf(self.check_limit_data(y_model)): 
                return -np.inf

        lh=-np.sum(((self.y_data-y_model)/self.y_err)[self.real_data_arg]**2)
        
        if np.isnan(lh): 
            return -np.inf
        return lh
    def run_mcmc(self):
        import emcee
        from multiprocessing import Pool
        from datetime import datetime
        import pickle
        np.random.seed(self.seed)
        pos =  self.p0+ 1e-3 * np.random.random((self.nwalker,self.ndim))
        # pos =  self.p0+ 1e-3 * np.random.randn(self.nwalker,self.ndim)

        if self.save_chain:
            # 将时间格式化为字符串，精确到分钟
            time_str = datetime.now().strftime('%Y-%m-%d %H:%M').replace(" ","")
            self.output_name=filename=f"{self.func_name}_emcee_output_{time_str}.h5"
            if os.path.exists(filename):
                os.remove(filename)
            backend = emcee.backends.HDFBackend(filename)
        
            with Pool(self.N_pool) as pool:
                sampler = emcee.EnsembleSampler(self.nwalker, 
                                                    self.ndim, 
                                                    self.log_prob, 
                                                    args=(self.args),
                                                    pool=pool,
                                                    backend=backend)
                _=sampler.run_mcmc(pos,self.nsamplers, progress=True)
                print(self.output_name,"has been saved in",os.path.abspath(os.path.curdir))


        else:
            with Pool(self.N_pool) as pool:
                sampler = emcee.EnsembleSampler(self.nwalker, 
                                                    self.ndim, 
                                                    self.log_prob, 
                                                    args=(self.args),pool=pool)
                _=sampler.run_mcmc(pos,self.nsamplers, progress=True)

        self.samplers=sampler

    def display_latex_or_text(self,string):
        from IPython.display import display, Math
        try:
            display(Math(string))
        except:
            print(string)


    def display_h5results(h5name=None,ParamterLabel=[],discard=1000,corner_range=0.99):
        import emcee
        from corner import corner
        def display_latex_or_text(string):
            from IPython.display import display, Math
            try:
                display(Math(string))
            except:
                print(string)

        reader = emcee.backends.HDFBackend(h5name)
        flat_samples = reader.get_chain(discard=discard, thin=15,flat=True)
        ParamterNum  = len(flat_samples[0])
        FitResult    = np.array([np.percentile(flat_samples[:,i],[16,50,84]) for i in range(ParamterNum)])
        FitResultVal = FitResult[:,1]
        FitResultStd = np.diff(FitResult)
        if not any(list(ParamterLabel)):
            ParamterLabel = [r"$Param_{%d}$"%i for i in range(ParamterNum)]

        corner(flat_samples,
            labels=ParamterLabel,
            truths=FitResultVal,
            truth_color="k",quantiles=[0.16, 0.5, 0.84], 
            title_fmt='.2f',show_titles=True,
            range=[corner_range]*ParamterNum)
        
        for i,val in enumerate(FitResultVal):
            display_latex_or_text(f"${ParamterLabel[i]}={val:.2f} + {FitResultStd[i][0]:.2f} - {FitResultStd[i][1]:.2f}$")
        

    def display_results(self,ParamterLabel=[],corner_range=0.99):
        self.flat_samples = self.samplers.get_chain(discard=self.nburn, thin=15, flat=True)
        self.ParamterNum  = len(self.flat_samples[0])
        self.FitResult    = np.array([np.percentile(self.flat_samples[:,i],[16,50,84]) for i in range(self.ParamterNum)])
        self.FitResultVal = self.FitResult[:,1]
        self.FitResultStd = np.diff(self.FitResult)
        if not any(list(ParamterLabel)):
            import inspect
            parameters = inspect.signature(self.func).parameters
            ParamterLabel = [param.name for param in parameters.values()][1:]
            if len(ParamterLabel)!=self.ParamterNum:
                ParamterLabel = [r"$Param_{%d}$"%i for i in range(self.ParamterNum)]

        self.ParamterLabel=ParamterLabel
        from corner import corner
        corner(self.flat_samples,
            labels=ParamterLabel,
            truths=self.FitResultVal,
            truth_color="k",quantiles=[0.16, 0.5, 0.84], 
            title_fmt='.2f',show_titles=True,
            range=[corner_range]*self.ParamterNum)
        
        for i,val in enumerate(self.FitResultVal):
            l=ParamterLabel[i].strip("$")
            self.display_latex_or_text(f"${l}={val:.2f} _ {{-{self.FitResultStd[i][0]:.2f}}} ^ {{+{self.FitResultStd[i][1]:.2f}}}$")

    @property
    def best_model(self):
        return self.func(self.FitResultVal,self.x_data,self.args)
    @property
    def residual(self):
        return self.best_model -self.y_data
    @property
    def chi2(self):
        return np.sum((self.residual/self.y_err)**2)
    @property
    def reduced_chi2(self):
        return np.sum((self.residual/self.y_err)**2)/(len(self.x_data)-self.ndim)
    def latex_text(self,ParamterLabel=[]):
        if not any(ParamterLabel): ParamterLabel=self.ParamterLabel
        for i,val in enumerate(self.FitResultVal):
            print("%9s = $%.2f_{-%.2f}^{+%.2f}$"%(ParamterLabel[i],val,self.FitResultStd[i,0],self.FitResultStd[i,1]))


def CSP_II_photometry(SNtype,abspath=os.path.abspath("./"),):
    """
    SNtype:
    IIb Ib Ic Ic-BL 
    """
    import astropy.units as u
    CSP_II_photometry_path="/home/ysy/_data/source/2020adow/LC/Ibc_LC/data/Stritzinger_2018/CSP-I_SE_SN_photometry/ASCII_tables/"
    """
    Tab2-3.ascii.dat:
    standard star
    Tab4.ascii.dat:
    Optical photometry of SE SNe in the Swope natural system.

    Tab5.ascii.dat:
    Optical photometry of SE SNe in the du Pont natural system.

    Tab6.ascii.dat:
    NIR photometry of SE SNe in the Swope natural system

    Tab7.ascii.dat:
    NIR photometry of SE SNe in the Swope natural system.

    Tab8.ascii.dat:
    NIR photometry of SE SNe in the du Pont natural system.
    """

    """
    1-  6  A6    ---     SN        SN name
    8-  9  I2    h       RAh       Right ascension (J2000)
    11- 12  I2    min     RAm       Right ascension (J2000)
    14- 18  F5.2  s       RAs       Right ascension (J2000)
        20  A1    ---     DE-       Declination sign (J2000)
    21- 22  I2    deg     DEd       Declination (J2000)
    24- 25  I2    arcmin  DEm       Declination (J2000)
    27- 30  F4.1  arcsec  DEs       Declination (J2000)
    32- 54  A23   ---     Host      Host galaxy name
    56- 74  A19   ---     Morph     Morphology (1)
    76- 81  F6.4  ---     z         Heliocentric redshift
        82  A1    ---   n_z         [y] Note on z (2)
    84- 92  A9    ---     RDisc     Discovery reference
    94-111  A18   ---     Disc      Discovery group (3)
    113-118  A6    ---     Type      Spectral type
    120-122  I3    ---     Phase1    Phase range
    124-127  I4    ---     Phase2    Phase range
    129-135  F7.2  d       T(B)max   Time of maximum brightness
    137-140  F4.2  d     e_T(B)max   rms uncertainty on T(B)max
        141  A1    ---   n_T(B)max   [d] Note on T(B)max (4)
        143  I1    ---     T5        ? indicates optical photometry in du Pont
                                        natural system (in table5)
        145  I1    ---     T6        ? indicates NIR photometry in Swope natural
                                        system (in table6)
        147  I1    ---     T7        ? indicates NIR photometry in Swope natural
                                        system (in table7)
        149  I1    ---     T8        ? indicates NIR photometry in du Pont natural
                                        system (in table8)
    """

    data_path_Swop  =os.path.join(abspath,CSP_II_photometry_path,"Tab4.ascii.dat")
    data_path_duPont=os.path.join(abspath,CSP_II_photometry_path,"Tab5.ascii.dat")

    SN_info_txt=open("/home/ysy/_data/source/2020adow/LC/Ibc_LC/data/Stritzinger_2018/SN_info.txt","r").readlines()[3:]
    SN_info_str=[i.strip("\n").split("\t") for i in SN_info_txt]
    SN_info_str=[[j.strip(" ")for j in i]for i in SN_info_str]
    SN_info_dict={i[0]:i[1:] for i in SN_info_str}
    
    
    data_Swop  =[i.strip().split() for i in open(data_path_Swop  ,"r").readlines()[3:]]
    data_duPont=[i.strip().split() for i in open(data_path_duPont,"r").readlines()[3:]]


    transient_name   =[i[0] for i in data_Swop   if len(i)<2]
    transient_name_dP=[i[0] for i in data_duPont if len(i)<2]
    
    transient_name=np.unique(transient_name)
    
    transient_index_Swop  =np.array([data_Swop  .index([i]) for i in transient_name if [i] in data_Swop])
    transient_index_duPont=np.array([data_duPont.index([i]) for i in transient_name if [i] in data_duPont])
    
    transient_index_duPont=transient_index_duPont[np.argsort(transient_index_duPont)]
    transient_index_Swop  =transient_index_Swop  [np.argsort(transient_index_Swop)]
    
    transient_index_Swop  =np.append(transient_index_Swop,-1)
    transient_index_duPont=np.append(transient_index_duPont,-1)

    CSP_ii_LC={}
    allband="ugriBV"
    for i in range(len(transient_name)):
        n=transient_name[i]
        
        if SNtype in SN_info_dict[n][2]:
        # print(n)
            distmod=Distance(float(SN_info_dict[transient_name[i]][10].split()[0])*u.Mpc).distmod.value
            Rv_host=float(SN_info_dict[n][7]) if SN_info_dict[n][7]!="..." else 3.1
            Av_MW  =float(SN_info_dict[n][6]) if SN_info_dict[n][6]!="..." else 0.0
            Av_host    =float(SN_info_dict[n][8].split("(")[0])               if (SN_info_dict[n][8]!="..." and SN_info_dict[n][8]!="0") else 0.0
            Av_host_err=float(SN_info_dict[n][8].split("(")[1].split(")")[0]) if (SN_info_dict[n][8]!="..." and SN_info_dict[n][8]!="0") else 0.0

            ebv_MW = Av_MW/3.1

            ebv_host    =Av_host/Rv_host
            ebv_host_err=Av_host_err/Rv_host
            # ebv_MW  =float(SN_info_dict[n][8].split("(")[0]) if (SN_info_dict[n][8]!="..." and SN_info_dict[n][8]!="0") else 0.0
            CSP_ii_LC[n]=LightCurve(name=n,
                                    band=allband,
                                    distmod=distmod,
                                    ebv_MW  =[ebv_MW,0.00],
                                    # Rv_host =Rv_host,
                                    ebv_host=[ebv_host,ebv_host_err],
                                    )
            
            data=data_Swop[transient_index_Swop[i]+1:transient_index_Swop[i+1]]
            data=np.array([[ float(d) if (d!="INDEF" and d!="UBDEF") else np.nan for d in D ]for D in data])
            data_resort=[]
            for d in data:
                for b in range(len(allband)):
                    data_resort.append([d[0],d[b*2+1],d[b*2+2],"Swop",allband[b],1])
            for  f in allband:
                CSP_ii_LC[n].__dict__[f].add_LightCurve_data(data_resort)
                CSP_ii_LC[n].__dict__[f]._mjd=CSP_ii_LC[n].__dict__[f]._mjd-2400000.5
    # for i in range(len(transient_name_dP)):
    #     n=transient_name[i]
    #     data=data_duPont[transient_index_duPont[i]+1:transient_index_duPont[i+1]]
    #     data=np.array([[ float(d) if (d!="INDEF" and d!="UBDEF") else np.nan for d in D ]for D in data])
    #     data_resort=[]
    #     for d in data:
    #         for b in range(len(allband)):
    #             data_resort.append([d[0],d[b*2+1],d[b*2+2],"duPont",allband[b],1])
    #     for  f in allband:
    #         CSP_ii_LC[n].__dict__[f].add_LightCurve_data(data_resort)
    return CSP_ii_LC

def match_same_date(l1,l2):

    if type(l1)!=np.array:l1=np.array(l1)
    if type(l2)!=np.array:l2=np.array(l2)

    min_gap=0.3
    l1_1night_all=[]
    for i1 in range(len(l1)):
        l1_1night=np.where(np.abs(l1-l1[i1])<min_gap)[0]
        l1_1night_all.append(list(l1_1night))
    
    All_l1_1night_index=[]
    [All_l1_1night_index.append(i) for i in l1_1night_all if i not in All_l1_1night_index]

    
    match_index=[]
    for i1 in range(len(All_l1_1night_index)):
        mean_time_l1=np.mean(l1[np.array(All_l1_1night_index[i1])])
        l2_1night_index=list(np.where(np.abs(l2-mean_time_l1)<min_gap)[0])
        if any(l2_1night_index):
            match_index.append([All_l1_1night_index[i1],l2_1night_index])
    return match_index
    

def colorLC(lc1,lc2):
    try:
        lc1=np.array(list(zip(lc1.mjd,lc1.mag,lc1.mag_e)))[lc1.detection==1,:]
        lc2=np.array(list(zip(lc2.mjd,lc2.mag,lc2.mag_e)))[lc2.detection==1,:]
    except:
        lc1=np.array(list(zip(lc1.mjd,lc1.mag,lc1.mag_e)))
        lc2=np.array(list(zip(lc2.mjd,lc2.mag,lc2.mag_e)))
    match_index=match_same_date(lc1[:,0],lc2[:,0])

    mjd1=np.array([np.mean(lc1[:,0][np.array(i[0])]) for i in match_index])
    mjd2=np.array([np.mean(lc2[:,0][np.array(i[1])]) for i in match_index])

    mag1=np.array([np.mean(lc1[:,1][np.array(i[0])]) for i in match_index])
    mag2=np.array([np.mean(lc2[:,1][np.array(i[1])]) for i in match_index])
    
    mage1=np.array([np.mean((lc1[:,2][np.array(i[0])])**2)**0.5 for i in match_index])
    mage2=np.array([np.mean((lc2[:,2][np.array(i[1])])**2)**0.5 for i in match_index])
    
    mjd=(mjd1+mjd2)/2
    color=mag1-mag2
    colore=(mage1**2+mage2**2)**0.5
    
    argsort=np.argsort(mjd)
    mjd    =mjd[argsort]
    color  =color[argsort]
    colore =colore[argsort]
    # names=["mjd","color","colorErr"]
    # formats=["float_"]*3
    # ],dtype=np.dtype({'names': names, 'formats': formats})

    array_mask=~np.isnan(color)*~np.isnan(colore)
    mjd    =mjd  [array_mask]
    color  =color[array_mask]
    colore =colore[array_mask]
    result=np.array([[i[0],i[1],i[2],] for i in zip(mjd,color,colore)])
    return result
class MCMC_Vacca_Leibundgut():
    def __init__(self,unit="mag"):
        # self.order=len(init_theta)-1
        # self.init_theta=init_theta
        # self.order=5
        # self.range={}
        # self.MIN_log_Prb=-1000
        # self.all_param_range=[[-np.inf,np.inf]]*(self.order+1)
        self.unit  =unit
        self.punish_UL=1000.0

    def calculate_all_LC_peak_mjd(self,LightCurve,sampler=1000,fraction=0.8,Boostrap=False,seed=100):
        LC_property_all={}
        peak_date_all={}
        for fitband in LightCurve.allband:
            LC_property=self.LC_property(LightCurve.__dict__[fitband],
                                            sampler=sampler,
                                            fraction=fraction,
                                            Boostrap=Boostrap,
                                            seed=seed)
            # peak_date,peak_mag,x_model,y_model,is_maximum,fit_parameter=self.run_curvefit(
            #                             x   =lc.__dict__[fitband].mjd,
            #                             y   =lc.__dict__[fitband].mag,
            #                             yerr=lc.__dict__[fitband].mag_e,
            #                             ydet=lc.__dict__[fitband].detection,)

            
            
            LightCurve.__dict__[fitband].peak_mjd = LC_property["peak_date"]
            LightCurve.__dict__[fitband].Vacca_Leibundgut = LC_property
            # lc.__dict__[fitband].Vacca_Leibundgut["LC_time"    ]=x_model
            # lc.__dict__[fitband].Vacca_Leibundgut["LC_mag"     ]=y_model
            # lc.__dict__[fitband].Vacca_Leibundgut["peak_mag"   ]=peak_mag-lc.distmod
            # lc.__dict__[fitband].Vacca_Leibundgut["param"      ]=fit_parameter
            # lc.__dict__[fitband].Vacca_Leibundgut["is_peak"    ]=is_maximum
            # lc.__dict__[fitband].Vacca_Leibundgut["decay_time" ]=decay_time
            # lc.__dict__[fitband].Vacca_Leibundgut["delta_mag15"]=delta_mag15
            

            peak_date_all[fitband]=(LC_property["peak_date"],
                                    LC_property["peak_mag"],
                                    LC_property["X_model"],
                                    LC_property["Y_model"],
                                    LC_property["is_maximum"],
                                    LC_property["fit_parameter_all"])
            LC_property_all[fitband]=LC_property
        return LC_property_all
        # return peak_date_all
        # return



    def Vacca_Leibundgut_LC(self,t,texp,t0,y0,m,g0,sigma0,theta,tao):
        # =tuple(theta)
        # texp,t0,y0,m,theta,tao=tuple(theta)
        # texp=theta[0]

        t = t-texp
        Mag  = y0+m*(t-t0)
        Mag += g0*np.exp(-(t-t0)**2/(2*sigma0**2))
        Mag /= 1-np.exp((tao-t)/theta)
        # Mag  = (y0+m*(t-t0)+g0*np.exp(-(t-t0)**2/(2*sigma0**2)))/(1-np.exp((tao-t)/theta))
        # plt.plot(t,10*(1-np.exp((tao-t)/theta)))
        # plt.plot(t,(g0*np.exp(-(t-t0)**2/(2*sigma0**2))))
        # plt.plot(t,y0+m*(t-t0))
        return Mag

    
    def log_LikeliHood(self,theta,x,y,yerr):
        x=np.array(x)
        y=np.array(y)
        yerr=np.array(yerr)

        y_model= self.Vacca_Leibundgut_LC(t,theta)
        sigma2 = yerr ** 2
        LikelyHood = -0.5 * np.sum((y - y_model) ** 2 / sigma2 + np.log(sigma2))
        return LikelyHood

    def Range(self):
    # TODO 
    # def Range(self,*arg):
        # arg=((1,2),(2,3))

        self.all_param_range={
            # "texp":[0,np.inf],
            "texp":[-20,0],
            "t0":[0,40],
            "y0":[-30,30],
            "m" :[0,30],
            "g0":[-10,0],
            "sigma0":[2,40],
            "theta":[-10,30],
            "tao":[-50,-1],
            # [0,10,17,0.05,-2,8,5,-10,]
        }
        self.all_param_range=np.array(list(self.all_param_range.values())).T
    # def log_prior(theta,*arg):


    
    def log_prior(self,theta):
        In_Range=True
        for r in range(len(theta)):
            # print(self.all_param_range[r][0],theta[r],self.all_param_range[r][1])
            if self.all_param_range[r][0]>theta[r] or \
               self.all_param_range[r][1]<theta[r]:
                In_Range=False
                # print(True)
                # print(In_Range)
        if theta[1]*theta[3]+theta[2]>0:
            pass
        else:
            In_Range=False

        if In_Range:
            return 0.0
        else:
            return -np.inf

    def log_uplimit(self,theta,x,y):
        if self.unit not in ["flux","mag"]:
            print("unit false self.unit should be 'mag' or 'flux'")
        Whether={"flux":"","mag":"not "}
        chi2=0
        for i in range(len(x)):
            
            y_modle=self.Vacca_Leibundgut_LC(theta,x[i])
            if y_modle > y[i]:
                model_beyond_data=Whether[self.unit]+"True"
            else:
                model_beyond_data=Whether[self.unit]+"False"
            
            model_beyond_data_bool=eval(model_beyond_data)

            if model_beyond_data_bool:
                chi2 = chi2 + (y[i] / y_modle * self.punish_UL)**2

        return -chi2
    
    def log_probability(self,theta,x,y,yerr,ydet):
        non_det=np.where(ydet==0)
        det    =np.where(ydet==1)
        #分类

        Xup    =x[non_det]
        Yup    =y[non_det]
        Yreal  =   y[det]
        Xreal  =   x[det]
        Yreale =yerr[det]
        log_UpL=self.log_uplimit(theta,Xup,Yup)
        log_Pr =self.log_prior  (theta)
        log_PrB=log_Pr+log_UpL+self.log_LikeliHood(theta, Xreal, Yreal, Yreale)
        return log_PrB
    
    def run_curvefit(self,x,y,yerr,ydet,initial_guess=[0,10,17,0.05,-2,8,5,-10,]):
        import matplotlib.pyplot as plt
        from scipy.optimize import curve_fit

        if type(x   )!=np.array:x=np.array(x)
        if type(y   )!=np.array:y=np.array(y)
        if type(yerr)!=np.array:yerr=np.array(yerr)
        if type(ydet)!=np.array:ydet=np.array(ydet)

        if not any(x):
            return np.nan,np.nan,np.array([]),np.array([]),False,np.array([])
        
        mjdindex=x.argsort()
        x=x[mjdindex]
        y=y[mjdindex]
        yerr=yerr[mjdindex]
        ydet=ydet[mjdindex]
        # print(y)
        mjdindex=np.where(~np.isnan(y))
        x=x[mjdindex]
        y=y[mjdindex]
        yerr=yerr[mjdindex]
        ydet=ydet[mjdindex]

        det_arg=np.where(ydet==1)
        x=x[det_arg]
        y=y[det_arg]
        yerr=yerr[det_arg]
        xshift=np.nanmin(x)
        
        # initial_guess[0]=np.nanmin(self.all_param_range[:,0])
        initial_guess[0]=self.all_param_range[1,0]
        yshift=16-np.nanmin(y) if any(y<0) else 0
        x=x-xshift
        # print(mjdindex)
        # print(x,y)
        if len(y)>=len(self.all_param_range):
            
            FitResult,FitCov=curve_fit( f=self.Vacca_Leibundgut_LC,
                                        xdata=x,
                                        ydata=y+yshift,
                                        sigma=yerr,
                                        p0=initial_guess,
                                        maxfev = 1800000,
                                        bounds=self.all_param_range,
                                        )#
            
            # print(FitResult)
            # print()
            # x_model_0=np.linspace(FitResult[0],np.nanmax(x),1000)
            x_model_0=np.linspace(FitResult[0],np.nanmax(x),1000)
            y_model_0=self.Vacca_Leibundgut_LC(x_model_0,*tuple(FitResult))-yshift
            peak_mag =np.nanmin(y_model_0)
            peak_date=x_model_0[np.where(y_model_0==peak_mag)]+xshift
            peak_date=peak_date[0]
            if peak_date==x_model_0[0]+xshift or peak_date==x_model_0[-1]+xshift:
                return peak_date,peak_mag,x_model_0+xshift,y_model_0,False,FitResult
            else:
                # x_model=np.linspace(np.nanmax([FitResult[0],-0.02*np.nanmax(x)]),1.02*np.nanmax(x),1000)
                # y_model=self.Vacca_Leibundgut_LC(x_model,*tuple(FitResult))
                # peak_date=x_model[np.where(y_model==np.nanmin(y_model))]+xshift
                return peak_date,peak_mag,x_model_0+xshift,y_model_0,True,FitResult
        else:
            return np.nan,np.nan,np.array([]),np.array([]),False,np.array([])
        
    # def run_mcmc(self,x,y,yerr,ydet,initial_guess=[]):
    #     import matplotlib.pyplot as plt
    #     mjdindex=x.argsort()
    #     x=x[mjdindex]
    #     y=y[mjdindex]
    #     yerr=yerr[mjdindex]
    #     ydet=ydet[mjdindex]
        
    #     import emcee
    #     # pos =  soln.x+ 1e-4 * np.random.randn(32, 3)
    #     np.random.seed(42)
    #     pos =  np.array(initial_guess)+ 1e-3 * np.random.randn(32,len(self.all_param_range))
    #     nwalkers, ndim = pos.shape

    #     import warnings
    #     #suppress warnings
    #     warnings.filterwarnings('ignore')        

    #     sampler = emcee.EnsembleSampler(nwalkers, ndim, self.log_probability, args=(x, y, yerr,ydet))
    #     _=sampler.run_mcmc(pos, 10000, progress=True)

    #     plt.figure(figsize=(15,8))
    #     plt.errorbar(x,y,yerr=yerr,fmt='.k')
    #     # plot_LC(lc,Xshift=max_mjd,Yshift=0)
    #     # plot_LC(LC.B,Xshift=max_mjd,Yshift=0)

    #     plt.legend(fontsize=14)
    #     plt.gca().invert_yaxis() 
    #     plt.xlabel("Date - $r_{max}$")
    #     plt.ylabel("y")
        
    #     flat_samples = sampler.get_chain(discard=4000, thin=15, flat=True)

    #     logPrb_sample=np.array([self.log_probability(t, x, y, yerr,ydet) for t in flat_samples])
    #     min_logPrb=-100
    #     min_logPrb=-0
    #     logPrb_sample_mask=np.where(logPrb_sample > min_logPrb)
    #     flat_samples_filtered= flat_samples[logPrb_sample_mask]

    #     x_resample=x
    #     x_resample=np.linspace(min(x),max(x),100)
    #     ReSample=[[self.Vacca_Leibundgut_LC(theta,xx) for theta in flat_samples_filtered]for xx in x_resample]
    #     ReSample_var  =np.array([np.std (i) for i in ReSample])
    #     ReSample_mean =np.array([np.mean(i) for i in ReSample])

    #     plt.fill_between(x_resample,ReSample_mean-ReSample_var , ReSample_mean+ReSample_var,
    #                 color="k", alpha=0.2)
    #     plt.plot(x_resample,ReSample_mean,"k")
        
    #     # plt.ylim(20,15)
    #     # return np.array(list(zip(ReSample_mean,ReSample_var))),flat_samples_filtered
    #     return flat_samples_filtered
    # def pred_curve(self,x,flat_samples):
    #     flat_samples= flat_samples
    #     x_resample=x
    #     ReSample=[[self.Vacca_Leibundgut_LC(theta,xx) for theta in flat_samples]for xx in x_resample]
    #     ReSample_var  =np.array([np.std (i) for i in ReSample])
    #     ReSample_mean =np.array([np.mean(i) for i in ReSample])
    #     return ReSample_mean,ReSample_var
    
    # def decay_rate(self,lc):
    #     lc_decay_rate={}
    #     for fitband in lc.allband:
    #         peak_mag=lc.__dict__[fitband].Vacca_Leibundgut["peak_mag"]
    #         half_peak_mag=peak_mag*0.75

    #         x=lc.__dict__[fitband].Vacca_Leibundgut["LC_time"]
    #         y=lc.__dict__[fitband].Vacca_Leibundgut["LC_mag" ]
    #         I=x>lc.__dict__[fitband].peak_mjd
    #         i=np.max(np.where(y[I]<peak_mag+0.75)[0])
    #         half_peak_mjd=x[I][i]
    #         decay_time=np.abs(half_peak_mjd-lc.__dict__[fitband].peak_mjd)
    #         lc_decay_rate[fitband]=decay_time
    #     return lc_decay_rate
    # def Bootstrap(self,LightCurve,sampler=1000,fraction=0.5,Boostrap=False):
    @property
    def LC_properties(self):
        return {'fit_parameter_all':np.array([[]]), 
                'is_maximum':False, 
                'X_model':np.array([]), 
                'Y_model':np.array([]), 
                'Y_model_err':np.array([]), 
                'peak_mag':np.nan, 
                'peak_mag_err':np.nan, 
                'peak_date':np.nan, 
                'peak_date_err_l':np.nan, 
                'peak_date_err_u':np.nan, 
                'rise_time':np.nan, 
                'rise_time_err_u':np.nan, 
                'rise_time_err_l':np.nan, 
                'delta_m15_rise':np.nan, 
                'delta_m15_err_rise':np.nan, 
                'half_peak_mjd_rise':np.nan, 
                'half_peak_mjd_rise_err_u':np.nan, 
                'half_peak_mjd_rise_err_l':np.nan, 
                'm15_rise_err':np.nan, 
                'decay_time':np.nan, 
                'decay_time_err_u':np.nan, 
                'decay_time_err_l':np.nan, 
                'delta_m15_decay':np.nan, 
                'delta_m15_err_decay':np.nan, 
                'half_peak_mjd_decay':np.nan, 
                'half_peak_mjd_decay_err_u':np.nan, 
                'half_peak_mjd_decay_err_l':np.nan, 
                'm15_decay_err':np.nan, 
                'm15':np.nan}
    def LC_property(self,LightCurve,sampler=1000,fraction=0.5,Boostrap=False,seed=100):

        LC_properties=self.LC_properties
        if not any(LightCurve.mjd):
            return LC_properties

        import warnings 
        warnings.filterwarnings("ignore")
        from tqdm import tqdm
        
        
        print("band:",LightCurve.filter if "filter" in LightCurve.__dict__ else "")
        
        _,_,x_model_0,y_model_0,is_maximum_0,fit_parameter_0=self.run_curvefit(
                                                x   =LightCurve.mjd,
                                                y   =LightCurve.mag,
                                                yerr=LightCurve.mag_e,
                                                ydet=LightCurve.detection,)
        
        # y_model_all.append(y_model)
        # x_model_all.append(x_model)
        # fit_parameter_all.append(fit_parameter)
        # fit_parameter_all.append(fit_parameter)


        y_model_BS=[]
        x_model_BS=[]
        fit_parameter_BS=[]

        y_model_BS.append(y_model_0)
        x_model_BS.append(x_model_0)
        fit_parameter_BS.append(fit_parameter_0)



        if Boostrap:
            min_texp_mjd=np.nanmin(LightCurve.mjd[LightCurve.detection==1])
            # print()
            N_fit_data=int(len(LightCurve.mjd)*fraction) if len(LightCurve.mjd)>30 else len(LightCurve.mjd)-1
            # N_fit_data=len(LightCurve.mjd)-1
            np.random.seed(seed)
            with tqdm(total=sampler, desc="Processing", unit="iteration") as pbar:
                for i in list(range(sampler)):
                # for i in list(range(sampler)):
                    
                    I=np.random.randint(0,len(LightCurve.mjd),N_fit_data)
                    self.all_param_range[1,0]=min_texp_mjd-np.nanmin(LightCurve.mjd[I][(LightCurve.detection[I]==1)*\
                                                                                       (~np.isnan(LightCurve.mag[I]))])
                    # print(self.all_param_range[:,0])
                    _,_,x_model,y_model,is_maximum,fit_parameter=self.run_curvefit(
                                                x   =LightCurve.mjd[I],
                                                y   =LightCurve.mag[I],
                                                yerr=LightCurve.mag_e[I],
                                                ydet=LightCurve.detection[I],)
                    
                    y_model_BS.append(y_model)
                    x_model_BS.append(x_model)
                    fit_parameter_BS.append(fit_parameter)
                    pbar.update(1)
        # else:
            
        self.all_param_range[1,0]=0

        fit_parameter_all=np.array(fit_parameter_BS)
        X_model=np.linspace(np.nanmin(fit_parameter_0[0])+np.nanmin(LightCurve.mjd),
                            np.nanmax(LightCurve.mjd),10000)
        Y_model_all=[]
        for i in range(len(y_model_BS)):
            y_temp=np.interp(X_model,x_model_BS[i],y_model_BS[i],left=np.nan)
            y_temp[np.where(np.isnan(y_temp))]=np.nan
            y_temp[np.where((y_temp>100) | (y_temp<-100))]=np.nan
            Y_model_all.append(y_temp)
                
        # y_model_mean=np.nanmean(Y_model_all,axis=0)
        Y_model_all =np.array(Y_model_all)        
        y_model_mean=np.interp(X_model,x_model_0,y_model_0,left=np.nan)

        def calculate_errors(y0_mean, y_data):

            # 计算y0的标准差
            y0_std_above = np.zeros_like(y0_mean)
            y0_std_below = np.zeros_like(y0_mean)
            for i in range(len(y0_mean)):
                # 找到每个x对应的y值
                y_values = y_data[:, i]
                # 计算均值之上的数据标准差
                y0_std_above[i] = 2*np.nanstd(y_values[y_values > y0_mean[i]])
                # 计算均值之下的数据标准差
                y0_std_below[i] = 2*np.nanstd(y_values[y_values < y0_mean[i]])
            return y0_std_above, y0_std_below
        y_model_std_l,y_model_std_u=calculate_errors(y_model_mean, Y_model_all[1:])
        I_argsort   =np.argsort(X_model)
        X_model     =X_model[I_argsort]
        y_model_mean=y_model_mean[I_argsort]
        y_model_std_l =y_model_std_l [I_argsort]
        y_model_std_u =y_model_std_u [I_argsort]
        

        peak_mag    =np.nanmin(y_model_mean)
        index_peak  =np.where(y_model_mean==peak_mag)[0][0]
        peak_mag_err_l=y_model_std_l[index_peak]
        peak_mag_err_u=y_model_std_u[index_peak]
        peak_date   =X_model    [index_peak]

        if np.where(X_model<peak_date)[0].any() and np.where(X_model>peak_date)[0].any():
            if np.where(np.isnan(X_model[X_model<peak_date]))[0].shape[0]==len(X_model[X_model<peak_date]) or \
               np.where(np.isnan(X_model[X_model>peak_date]))[0].shape[0]==len(X_model[X_model>peak_date]):
                is_maximum=False
            else:
                is_maximum=True

        else:
            is_maximum=False

        
        if not is_maximum:
            LC_properties['peak_mag']=peak_mag
            LC_properties['peak_mag_err_l']=peak_mag_err_l
            LC_properties['peak_mag_err_u']=peak_mag_err_u
            LC_properties['peak_date']=peak_date
            LC_properties['is_maximum']=is_maximum
            LC_properties['X_model']=X_model
            LC_properties['Y_model']=y_model_mean
            LC_properties['Y_model_err_l']=y_model_std_l
            LC_properties['Y_model_err_u']=y_model_std_u
            LC_properties["fit_parameter_all"]=fit_parameter_all
            return LC_properties

        # print(X_model[:10])
        # print(peak_date)
        # print(X_model[(y_model_mean-y_model_std)<peak_mag])
        # print(y_model_mean-y_model_std)
        # print(y_model_std)
        # print(y_model_mean)
        def split_array(index,index_peak):
            # 在采样的数据中，早期光变的起伏非常大，所以需要用以下这个function去排除掉早期的光变采样
            result = []
            sublist = [index[0]]

            for i in range(1, len(index)):
                if index[i] - sublist[-1] > 1:
                    result.append(sublist)
                    sublist = [index[i]]
                else:
                    sublist.append(index[i])

            result.append(sublist)
            for i in range(len(result)):
                if index_peak in result[i]:
                    return np.array(result[i])
            return np.array(index)
        
        i_PrePeak=np.where((y_model_mean-y_model_std_u)<peak_mag)[0]
        i_PrePeak=split_array(i_PrePeak,index_peak) if Boostrap else np.array([])
        peak_date_err_l =np.abs(np.nanmin(X_model[i_PrePeak])-peak_date) if Boostrap else 0
        peak_date_err_u =np.abs(np.nanmax(X_model[i_PrePeak])-peak_date) if Boostrap else 0

        
        print("peak date:%.2f -%.2f +%.2f"%(peak_date,peak_date_err_l,peak_date_err_u))
        print("peak mag:%.2f -%.2f +%.2f"%(peak_mag,peak_mag_err_l,peak_mag_err_u))
      
        LC_properties["fit_parameter_all"]=fit_parameter_all
        LC_properties["is_maximum"] =True if Boostrap else is_maximum
        LC_properties["X_model"]    =X_model
        LC_properties["Y_model"]    =y_model_mean
        LC_properties["Y_model_err_u"]=y_model_std_u
        LC_properties["Y_model_err_l"]=y_model_std_l
        LC_properties["peak_mag"]   =peak_mag
        LC_properties["peak_mag_err_l"]=peak_mag_err_l
        LC_properties["peak_mag_err_u"]=peak_mag_err_u
        LC_properties["peak_date"]  =peak_date
        LC_properties["peak_date_err_l"]=peak_date_err_l
        LC_properties["peak_date_err_u"]=peak_date_err_u


        # pre peak
        I  =X_model<(peak_date)
        I15=X_model>(peak_date-15)

        m15    =y_model_mean[I*I15][0]
        m15_err_l=y_model_std_l[I*I15][0]
        m15_err_u=y_model_std_u[I*I15][0]
        
        delta_m15    =np.abs(m15-peak_mag)
        delta_m15_err_u=np.nanmean(m15_err_l**2+peak_mag_err_u**2)**0.5
        delta_m15_err_l=np.nanmean(m15_err_u**2+peak_mag_err_l**2)**0.5
        
        i=np.nanmin(np.where(y_model_mean[I]<peak_mag+0.75)[0])
        i_PrePeak=np.where((y_model_mean[I]-y_model_std_u[I])<(peak_mag+0.75))[0]
        i_PrePeak=split_array(i_PrePeak,i)  if Boostrap else np.array([])

        half_peak_mjd= X_model[I][i]
        if Boostrap:
            half_peak_mjd_err_l =np.abs(X_model[I][i_PrePeak][0]-half_peak_mjd)   
            half_peak_mjd_err_u =np.abs(X_model[I][(y_model_mean[I]+y_model_std_l[I])<(peak_mag+0.75)][0]-half_peak_mjd) if peak_mag_err_l<0.75 else np.abs(half_peak_mjd-peak_date) 
        else:
            half_peak_mjd_err_l =np.nan
            half_peak_mjd_err_u =np.nan
                                
        

        rise_time       = np.abs(half_peak_mjd-peak_date)
        rise_time_err_u = np.nanmean(half_peak_mjd_err_l**2+peak_date_err_u**2)**0.5
        rise_time_err_l = np.nanmean(half_peak_mjd_err_u**2+peak_date_err_l**2)**0.5

        
        print("Rise time from half peak to peak:%.2f -%.2f +%.2f"%(rise_time,rise_time_err_u,rise_time_err_l))
        print("delta m15 (rise):%.2f -%.2f +%.2f"%(delta_m15,delta_m15_err_l,delta_m15_err_u))

        LC_properties["rise_time"]       =rise_time
        LC_properties["rise_time_err_u"] =rise_time_err_u
        LC_properties["rise_time_err_l"] =rise_time_err_l
        LC_properties["delta_m15_rise"]  =delta_m15
        LC_properties["delta_m15_rise_err_l"]=delta_m15_err_l
        LC_properties["delta_m15_rise_err_u"]=delta_m15_err_u

        LC_properties["half_peak_mjd_rise"]  =half_peak_mjd
        LC_properties["half_peak_mjd_rise_err_u"]  =half_peak_mjd_err_u
        LC_properties["half_peak_mjd_rise_err_l"]  =half_peak_mjd_err_l
        LC_properties["m15_rise_err_l"]  =m15_err_l
        LC_properties["m15_rise_err_u"]  =m15_err_u
        
        # post peak
        I=X_model>peak_date
        I15=X_model<(peak_date+15)


        m15    =y_model_mean[I*(X_model<peak_date+15)][-1]
        # m15_err=y_model_std [I*(X_model<peak_date+15)][-1]
        m15_err_l=y_model_std_l[I*I15][0]
        m15_err_u=y_model_std_u[I*I15][0]

        delta_m15    =np.abs(m15-peak_mag)
        # delta_m15_err=np.nanmean(m15_err**2+peak_mag_err**2)**0.5
        delta_m15_err_u=np.nanmean(m15_err_l**2+peak_mag_err_u**2)**0.5
        delta_m15_err_l=np.nanmean(m15_err_u**2+peak_mag_err_l**2)**0.5


        i=np.max(np.where(y_model_mean[I]<peak_mag+0.75)[0])
        half_peak_mjd= X_model[I][i]
        

        if Boostrap:
            half_peak_mjd_err_u =np.abs(X_model[I][(y_model_mean[I]-y_model_std_u[I])<peak_mag+0.75][-1]-half_peak_mjd) 
            half_peak_mjd_err_l =np.abs(X_model[I][(y_model_mean[I]+y_model_std_l[I])<peak_mag+0.75][-1]-half_peak_mjd) if peak_mag_err_l<0.75 else np.abs(half_peak_mjd-peak_date)           
        else:
            half_peak_mjd_err_l =np.nan
            half_peak_mjd_err_u =np.nan

        decay_time       = np.abs(half_peak_mjd-peak_date)
        decay_time_err_u = np.nanmean(half_peak_mjd_err_u**2+peak_date_err_l**2)**0.5
        decay_time_err_l = np.nanmean(half_peak_mjd_err_l**2+peak_date_err_u**2)**0.5
        
        print("Decay time from peak to half peak:%.2f -%.2f +%.2f"%(decay_time,decay_time_err_u,decay_time_err_l))
        print("delta m15(decay):%.2f -%.2f +%.2f"%(delta_m15,delta_m15_err_l,delta_m15_err_u))

        LC_properties["decay_time"]       =decay_time
        LC_properties["decay_time_err_u"] =decay_time_err_u
        LC_properties["decay_time_err_l"] =decay_time_err_l
        LC_properties["delta_m15_decay"]   =delta_m15
        LC_properties["delta_m15_decay_err_u"]=delta_m15_err_u
        LC_properties["delta_m15_decay_err_l"]=delta_m15_err_l
        LC_properties["half_peak_mjd_decay"]  =half_peak_mjd
        LC_properties["half_peak_mjd_decay_err_u"]  =half_peak_mjd_err_u
        LC_properties["half_peak_mjd_decay_err_l"]  =half_peak_mjd_err_l
        LC_properties["m15_decay_err_u"]  =m15_err_u
        LC_properties["m15_decay_err_l"]  =m15_err_l
        LC_properties["m15"]  =m15
        # print(LC_properties.keys())
        # print(y_model)        
        return LC_properties
    
class MagFlux():
    # 星等转流量 && 流量转星等
    import astropy.units as u
    # from SNeTools.instrument_info import zp_AB,zp_Vega,center_wv
    @property
    def AB(self):return "ab"
    @property
    def Vega(self):return "vega"
    @property
    def F_nu(self):return u.Jy
    @property
    def F_lambda(self):return u.Unit("erg cm-2 s-1 AA-1")
    # TODO add_filter_info
    
    def __init__(self,flux_unit=u.Unit("erg cm-2 s-1 AA-1"),default_system=True):
        import pandas as pd
        self.unit=flux_unit
        
        self.filter_info_csv_path="/home/ysy/_data/source/2020adow/model/filter_info.csv"

        self.filter_info_0=pd.read_csv(self.filter_info_csv_path,index_col=False)
        filter_info=pd.read_csv(self.filter_info_csv_path)

        filter_info["vega.flux.lambda"]=filter_info["vega.flux.lambda"]*1e-9
        filter_info["ab.flux.lambda"]  =filter_info["ab.flux.lambda"]*1e-9
        filter_info["ab.flux.Jy"]      =(filter_info[  "ab.flux.lambda"].values*self.F_lambda*\
                                         (filter_info["center_wv"].values*u.AA)**2/const.c).to_value(u.Jy)
        filter_info["vega.flux.Jy"]    =(filter_info["vega.flux.lambda"].values*self.F_lambda*\
                                         (filter_info["center_wv"].values*u.AA)**2/const.c).to_value(u.Jy)

        filter_info["ab.mag.lambda"]   =-2.5*np.log10(filter_info["ab.flux.lambda"])
        filter_info["ab.mag.Jy"]       =-2.5*np.log10(filter_info["ab.flux.Jy"])
        filter_info["vega.mag.Jy"]     =-2.5*np.log10(filter_info["vega.flux.Jy"])
        filter_info["vega.mag.lambda"] =-2.5*np.log10(filter_info["vega.flux.lambda"])

        filter_info["flux.lambda"]= [fi[fi["defualt_system"]+".flux.lambda"] for i,fi in filter_info.iterrows()] 
        filter_info["flux.Jy"   ] = [fi[fi["defualt_system"]+".flux.Jy"    ] for i,fi in filter_info.iterrows()]
        filter_info["mag.Jy"    ] = [fi[fi["defualt_system"]+".mag.Jy"    ] for i,fi in filter_info.iterrows()] 
        filter_info["mag.lambda"] = [fi[fi["defualt_system"]+".mag.lambda"] for i,fi in filter_info.iterrows()] 
        
        self.filter_info=filter_info
        self.defualt_system=default_system
    
        # magnitude_system=["ab","vega","default"]
        # if system.lower() not in magnitude_system:
        #     print("the specific magnitude system not in", magnitude_system)
        # else:
        #     self.system=system
    
      
    # def filters2_defualtsytem(self,filters):
        #########################################################
    def add_filter_info(self,new_info):
        all_columns=self.filter_info_0.columns
        if set(new_info)== set(all_columns):
            if new_info["filter"] in self.filter_info_0["filter"].values:
                import warnings
                warnings.warn(new_info["filter"]+" has already in filter_info")
                print(self.filter_info_csv_path)
                return self.filter_info_0
            else:
                new_info['defualt_system']=new_info['defualt_system'].lower()
                self.filter_info_0=self.filter_info_0.append(new_info,ignore_index=True)
                self.filter_info_0.to_csv(self.filter_info_csv_path,index=False)
                print("update in",self.filter_info_csv_path)
                return self.filter_info_0
        else:
            print("new_info should contain",all_columns)
            return
        
    def filters2center_wv(self,filters):
        if type(filters)==str or type(filters) == np.str_:
            if filters in self.filter_info_0["filter"].values:
                filters=[filters]
        for f in set(filters): 
            if f not in self.filter_info["filter"].values:
                from IPython.display import display
                display(self.filter_info_0)
                raise ValueError(f + "band is not in filter_info, please add. \n you can use self.add_filter_info() to add new info \n you can refer http://svo2.cab.inta-csic.es/theory/fps/")
        return self.array2float(np.array([self.filter_info[self.filter_info["filter"]==f]["center_wv"].values[0] for f in filters]) *u.AA)
    
    def filters2zp_mag_Jy(self,filters,photometry_system=[]):
        if type(filters)==str:
            if filters in self.filter_info["filter"].values:
                filters=[filters]
        if not any(photometry_system) or len(filters)!=len(photometry_system):
            return np.array([self.filter_info[self.filter_info["filter"]==f]["mag.Jy"].values[0] for f in filters])
        else:
            return np.array([self.filter_info[self.filter_info["filter"]==f][s+".mag.Jy"].values[0] for f,s in zip(filters,photometry_system)])

    def filters2zp_mag_lambda(self,filters,photometry_system=[]):
        if type(filters)==str:
            if filters in self.filter_info["filter"].values:
                filters=[filters]
        if not any(photometry_system) or len(filters)!=len(photometry_system):
            return np.array([self.filter_info[self.filter_info["filter"]==f]["mag.lambda"].values[0] for f in filters])
        else:
            return np.array([self.filter_info[self.filter_info["filter"]==f][s+".mag.lambda"].values[0] for f,s in zip(filters,photometry_system)])

    def filters2zp_flux_Jy(self,filters,photometry_system=[]):
        if type(filters)==str:
            if filters in self.filter_info["filter"].values:
                filters=[filters]
        if not any(photometry_system) or len(filters)!=len(photometry_system):
            return np.array([self.filter_info[self.filter_info["filter"]==f]["flux.Jy"].values[0] for f in filters])
        else:
            return np.array([self.filter_info[self.filter_info["filter"]==f][s+".flux.Jy"].values[0] for f,s in zip(filters,photometry_system)])
    def filters2zp_flux_lambda(self,filters,photometry_system=[]):
        if type(filters)==str:
            if filters in self.filter_info["filter"].values:
                filters=[filters]
        if not any(photometry_system) or len(filters)!=len(photometry_system):
            return np.array([self.filter_info[self.filter_info["filter"]==f]["flux.lambda"].values[0] for f in filters])
        else:
            return np.array([self.filter_info[self.filter_info["filter"]==f][s+".flux.lambda"].values[0] for f,s in zip(filters,photometry_system)])
    
    # def fluxvalue2unit(self,flux,filters=[]):
    #     if self.unit==self.F_lambda:
    #         return flux*self.unit
    #     else:
    #         c_WV=self.filters2center_wv(filters)*u.AA
    def array2float(self,array):
        if len(array)<=1:return array[0]
        return array
    def check_fluxunit(self,flux):
        if isinstance(flux,u.Quantity):
            if flux.unit==self.F_lambda or flux.unit==self.F_nu:
                return flux.unit
            else:
                from astropy.units import UnitTypeError
                raise UnitTypeError(flux.unit)
        else:
            from astropy.units import UnitsWarning
            UnitsWarning("flux have no units, and set as ",self.unit)            
            return self.unit
    def magnitude(self,magnitude,Filter,system=[]):
        if not (isinstance(magnitude, list) or isinstance(magnitude, np.ndarray)):
            magnitude=np.array([magnitude])
            Filter=np.array([Filter])
        if isinstance(Filter,np.str_) or isinstance(Filter,str):
            Filter=np.array([Filter])
        
        if self.unit==self.F_lambda: zp=self.filters2zp_flux_lambda(Filter,system)
        if self.unit==self.F_nu    : zp=self.filters2zp_flux_Jy    (Filter,system)
        flux= zp * 10**(magnitude/-2.5) * self.unit 
        
        return self.array2float(flux)
    
    def mage(self,magnitude,mage,Filter,system=[]):
        mage    = np.abs(mage)

        fluxe_u = self.magnitude(magnitude,Filter,system)-self.magnitude(magnitude+mage,Filter,system)
        fluxe   = fluxe_u
        return fluxe

    def flux(self,flux,Filter,system=[]):
        import astropy.constants as c
        import numpy as np


        unit=self.check_fluxunit(flux)
        if not (isinstance(Filter, list) or isinstance(Filter, np.ndarray)):
            Filter=np.array([Filter])
            flux  =np.array([flux.to_value(unit)])
        if isinstance(Filter,np.str_) or isinstance(Filter,str):
            Filter=np.array([Filter])
            
        if unit==self.F_lambda: zp=self.filters2zp_flux_lambda(Filter,system)
        if unit==self.F_nu    : zp=self.filters2zp_flux_Jy    (Filter,system)

        flux[flux<0]=np.nan
        if not isinstance(flux,np.ndarray):flux=flux.to_value()
        mag= -2.5 * np.log10(flux/zp)
        return self.array2float(mag)
    
class PlotSpec():
    def __init__(self, Spectra,max_spec=np.nan,wv_range=[2000,12000],figsize=[5,7],
                 ref_mjd=np.nan,
                 min_offset=0.3,max_offset=0.4,redshit=0,SNR_max=1):
        self.Spectra = Spectra
        self.max_spec = max_spec
        self.wv_range = wv_range
        self.figsize = figsize
        self.min_offset = min_offset
        self.max_offset = max_offset
        self.redshit = redshit
        self.SNR_max = SNR_max
        self.ref_mjd = ref_mjd if ~np.isnan(ref_mjd) else self.Spectra.ref_mjd
        self.spec_color_telescp = {
            "XLT": "b",
            "LJT" : "g",
            "ogg2m" : "r",
            "Keck"  : "purple",
            "Shane"  : "orange",
            "Ekar_AFOSC":"r",
            "Kast"  : "orange",
            "Unkown":"k",
            "NOT"   :"k",
            "LT"   :"coral",
            "HCT"  :"purple",
            "Other":"magenta",
            "P60"  :"olive", 
            "TNG":"darkred",
            "VLT":"skyblue",
        }
    def DER_SNR(self,flux):    
    # =====================================================================================
        """
        DESCRIPTION This function computes the signal to noise ratio DER_SNR following the
                    definition set forth by the Spectral Container Working Group of ST-ECF,
                MAST and CADC. 

                    signal = median(flux)      
                    noise  = 1.482602 / sqrt(6) median(abs(2 flux_i - flux_i-2 - flux_i+2))
                snr    = signal / noise
                    values with padded zeros are skipped

        USAGE       snr = DER_SNR(flux)
        PARAMETERS  none
        INPUT       flux (the computation is unit independent)
        OUTPUT      the estimated signal-to-noise ratio [dimensionless]
        USES        numpy      
        NOTES       The DER_SNR algorithm is an unbiased estimator describing the spectrum 
                as a whole as long as
                    * the noise is uncorrelated in wavelength bins spaced two pixels apart
                    * the noise is Normal distributed
                    * for large wavelength regions, the signal over the scale of 5 or
                    more pixels can be approximated by a straight line
        
                    For most spectra, these conditions are met.

        REFERENCES  * ST-ECF Newsletter, Issue #42:
                    www.spacetelescope.org/about/further_information/newsletters/html/newsletter_42.html
                    * Software:
                www.stecf.org/software/ASTROsoft/DER_SNR/
        AUTHOR      Felix Stoehr, ST-ECF
                    24.05.2007, fst, initial import
                    01.01.2007, fst, added more help text
                    28.04.2010, fst, return value is a float now instead of a numpy.float64
        """
        from numpy import array, where, median, abs 

        flux = array(flux)

        # Values that are exactly zero (padded) are skipped
        flux = array(flux[where(flux != 0.0)])
        n    = len(flux)      

        # For spectra shorter than this, no value can be returned
        if (n>4):
            signal = median(flux)

            noise  = 0.6052697 * median(abs(2.0 * flux[2:n-2] - flux[0:n-4] - flux[4:n]))

            return float(signal / noise)  

        else:

            return 0.0
    def spec_norm(self,spec,order=3):
        from scipy.optimize import curve_fit
        wv   = spec[:,0]
        flux = spec[:,1]
        def poly(x,*arg):
            y=np.poly1d(arg)
            return y(x)
        param,pcov=curve_fit(poly,xdata=wv,ydata=flux,p0=[1]*order)
        flux_ploy=poly(wv,*param)
        scale=2*np.nanmedian(flux_ploy)
        spec[:,1]=flux/scale
        return spec

    def spec_offset(self,spec1,spec2,min_offset=0.1,max_offset=0.3):
        # 波长覆盖不重叠
        if np.nanmax(spec1[:,0])<np.nanmin(spec2[:,0]) or \
           np.nanmin(spec1[:,0])>np.nanmax(spec2[:,0]) :
            return 0
        wv_range=np.array([ np.nanmin(spec1[:,0]),
                            np.nanmax(spec1[:,0]),
                            np.nanmin(spec2[:,0]),
                            np.nanmax(spec2[:,0]),
                ])
        # 波长重叠的部分
        wv_range  =wv_range[np.argsort(wv_range)][1:3]
        wv_intp   =np.linspace(wv_range[0],wv_range[1],len(spec1)+len(spec2))
        spec2_intp=np.interp(wv_intp,spec2[:,0],spec2[:,1])
        spec1_intp=np.interp(wv_intp,spec1[:,0],spec1[:,1])
        offset=np.nanmin(spec1_intp-spec2_intp)
        # print(offset)
        if offset<=0:
            return min_offset-offset
        elif offset>max_offset:
            return min_offset-offset
        elif min_offset<offset:
            return offset
        else:
            return 0
    def spec_cut(self,spec,wv_min,wv_max):
        spec_arg=np.where(((spec[:,0]>wv_min) &(spec[:,0]<wv_max) ))
        return spec[spec_arg]
    def spec_rebin(self,spec,step):
        # 将数据分成3个块，每个块包含6行
        block_size = step
        num_blocks = spec.shape[0] // block_size

        # 初始化一个3x7的结果数组
        result = np.zeros((num_blocks,spec.shape[1]))

        # 对每个块进行平均操作
        for i in range(num_blocks):
            start_row = i * block_size
            end_row = start_row + block_size
            block_data = spec[start_row:end_row, :]
            block_average = np.mean(block_data, axis=0)
            result[i, :] = block_average

        return result

    def spec_check_SNR(self,spec,SNR_max=20):
        useful=True
        spec_temp=spec.copy()
        k=2
        rebin=False
        while self.DER_SNR(spec[:,1])<SNR_max:
            spec=self.spec_rebin(spec_temp,k)
            rebin=k
            k+=1
            if k>15:
                useful=False
                break
        return spec,useful,rebin
    def plot_feature(center_wave,range=[-1000,1000],normalize_fit_range=[-2000,2000]):
        print()
    def plot_AllSpec_inst(self,lw=0.1,alpha=1,date_offset=100):
        name=self.Spectra.name
        spec_num = len(self.Spectra)

        offset=0
        spec_data2    =self.Spectra[0].wv_flux
        spec_data2    =self.spec_norm(spec_data2)
        spec_data2    =self.spec_cut(spec_data2,self.wv_range[0],self.wv_range[1])
        spec_data2,_,_=self.spec_check_SNR(spec_data2)
        # TODO spectra arrange by data
        # TODO smothh the spectra
        # TODO add the flux error

        
        
        plt.figure(figsize=(self.figsize[0],self.figsize[1]),dpi=150)
        # plot_color={}

        for i in range(spec_num)[:]:
            if ~np.isnan(self.max_spec):
                if (i%self.max_spec)==0:
                    plt.figure(figsize=(self.figsize[0],self.figsize[1]),dpi=100)
                
            spec_data1= spec_data2.copy()
            spec_data2= self.Spectra[i].wv_flux

            spec_data2= self.spec_norm(spec_data2)
            spec_data2= self.spec_cut (spec_data2,self.wv_range[0],self.wv_range[1])

            # 判断信噪比
            spec_data2,useful,rebin= self.spec_check_SNR(spec_data2,SNR_max=self.SNR_max)
            
            spec_mjd  = self.Spectra[i].time.mjd
            # if useful==False:
            #     print("%.2f(%d) is useless"%(spec_mjd-self.Spectra.ref_mjd,i))
            #     continue

            
            offset+=self.spec_offset(spec_data1,
                                     spec_data2,
                                     min_offset=self.min_offset,
                                     max_offset=self.max_offset)
            
            legend = plt.gca().get_legend()
            legend_labels=[]
            if legend:
                legend_labels = [text.get_text() for text in legend.get_texts()]
                
                # legend_labels = [text.get_text() for text in plt.gca().get_legend().get_texts()]
            if self.Spectra[i].telescope in legend_labels:
                plot_lable = None
            # plot_lable = "%.2f(%d)"%(spec_mjd-self.Spectra.ref_mjd,i)
            else:
                plot_lable = self.Spectra[i].telescope

            color=self.spec_color_telescp[self.Spectra[i].telescope]
            # plot_color[self.Spectra[i].telescope]=color

            if rebin!=False:LW=np.min([lw*rebin*0.2,1])
            else:LW=lw
            
            plt.plot(spec_data2[:,0]/(1+self.Spectra.redshift),
                     spec_data2[:,1]-offset,
                     label=plot_lable,
                     color=color,
                     lw=LW,
                     alpha=alpha,
                     )
            # if np.abs((spec_data1[-1,1]-offset)-(spec_data2[-1,1]-offset))<0.4:offset+0.4
            
            phase="%.1f"%(spec_mjd-self.Spectra.ref_mjd) if spec_mjd<self.ref_mjd else "+%.1f"%(spec_mjd-self.Spectra.ref_mjd)
            plt.text(spec_data2[-1,0]/(1+self.Spectra.redshift)+date_offset,
                     np.mean(spec_data2[-10:,1])-offset,
                     s=phase,
                     color=color,)
            
            
            # yshift_plot=1.3
            # xshift_plot=10000
            # k=0
            # for t in plot_color.keys():
            #     k+=0.2
            #     plt.text(200+xshift_plot, yshift_plot-k, t , fontsize=6, verticalalignment='center',color="k",alpha=1)
            #     plt.gca().add_patch(plt.Rectangle((-40+xshift_plot, yshift_plot-k-0.03), 200, 0.06, color=plot_color[t], alpha=1))
            
            
            legend = plt.gca().legend(bbox_to_anchor=(0.34, 0.92, 0.3, 0.2), loc='center',ncol=5,fontsize=8,frameon=False)
            for text in legend.get_texts():
                text.set_color('black')
            for handle in legend.legendHandles:
                handle.set_color(handle.get_color())
                handle.set_lw(8)
            # plt.xlim(self.wv_range[0]+date_offset,
            #          self.wv_range[1])
            # plt.yticks([])

    def plot_AllSpec(self,lw=0.1,alpha=1,date_offset=100):
        # name=self.Spectra.name
        spec_num = len(self.Spectra)

        offset=0
        spec_data2    =self.Spectra[0].wv_flux
        spec_data2    =self.spec_norm(spec_data2)
        spec_data2    =self.spec_cut(spec_data2,self.wv_range[0],self.wv_range[1])
        spec_data2,_,_=self.spec_check_SNR(spec_data2)
        # TODO spectra arrange by data
        # TODO smothh the spectra
        # TODO add the flux error

        plt.figure(figsize=(self.figsize[0],self.figsize[1]),dpi=150)

        for i in range(spec_num)[:]:
            if ~np.isnan(self.max_spec):
                if (i%self.max_spec)==0:
                    plt.figure(figsize=(self.figsize[0],self.figsize[1]),dpi=150)

                    plt.ylabel(r"Flux + offset ($erg~cm^{-2}~s^{-1}~\mathring{A}^{-1}$)")
                    plt.xlabel("Rest Wavelength (Angstrom)")
                
            spec_data1= spec_data2.copy()
            spec_data2= self.Spectra[i].wv_flux

            spec_data2= self.spec_norm(spec_data2)
            spec_data2= self.spec_cut (spec_data2,self.wv_range[0],self.wv_range[1])

            # 判断信噪比
            spec_data2,useful,rebin= self.spec_check_SNR(spec_data2,SNR_max=self.SNR_max)
            
            spec_mjd  = self.Spectra[i].time.mjd
            # if useful==False:
            #     print("%.2f(%d) is useless"%(spec_mjd-self.Spectra.ref_mjd,i))
            #     continue

            
            offset+=self.spec_offset(spec_data1,
                                     spec_data2,
                                     min_offset=self.min_offset,
                                     max_offset=self.max_offset)
            
            legend = plt.gca().get_legend()
                

            if rebin!=False:LW=np.min([lw*rebin*0.2,1])
            else:LW=lw
            phase="%.1f"%(spec_mjd-self.Spectra.ref_mjd) if spec_mjd<self.ref_mjd else "+%.1f"%(spec_mjd-self.Spectra.ref_mjd)            
            plt.plot(spec_data2[:,0]/(1+self.Spectra.redshift),
                     spec_data2[:,1]-offset,
                    #  label=plot_lable,
                    #  color=color,
                     lw=LW,
                     alpha=alpha,
                     label=f"{phase}({i})")
            
            legend = plt.gca().legend(bbox_to_anchor=(0.34, 0.98, 0.3, 0.2), loc='center',ncol=5,fontsize=8,frameon=False)
            for text in legend.get_texts():
                text.set_color('black')
            for handle in legend.legendHandles:
                handle.set_color(handle.get_color())
                handle.set_lw(8)
        

        
            # plt.xlim(self.wv_range[0]+date_offset,
            #          self.wv_range[1])
            # plt.yticks([])



class Spectra:
    def __init__(self,name,redshift=0.0,distmod=0,ebv_host=[np.nan,np.nan],ebv_MW=[np.nan,np.nan],ref_mjd=np.nan):
        self.name     =name
        self.redshift =redshift
        self.distmod  =distmod    
        self.ebv_host =ebv_host
        self.ebv_MW   =ebv_MW
        self.ref_mjd  =ref_mjd
        self.SpecAll  =[]
    def __len__(self):
        # 返回子对象的数量
        return len(self.SpecAll)
    def __getitem__(self,index):
        return self.SpecAll[index]
    # def generate_wavelength(self,ref_wavelength,ref_pixal,spec_resolution,spec_length):
    #     wv=np.linspace(\
    #     ref_wavelength-spec_resolution* (ref_pixal-1),\
    #     ref_wavelength-spec_resolution* (ref_pixal-spec_length),\
    #     spec_length)
    #     return wv
    # # def add_spec(self,spec_name=None):

    def readspec(self,spec_fits_name,time_key=None,time_format=None):
        from astropy.io import fits
        from astropy.wcs import WCS
        # 读取FITS文件

        from astropy.io import fits
        from astropy.wcs import WCS
        import numpy as np

        # 打开FITS文件
        spec_flux=fits.getdata(spec_fits_name)
        header=fits.getheader(spec_fits_name)

        # 从header中初始化WCS对象
        wcs = WCS(header)

        # 获取数据的形状
        # astropy 在读取wcs的时候，会读取头文件的时间。比如OBS-TIME OBS-MJD.这时候如果头文件格式不对，会造成很多warning，如何去掉warning
        shape = spec_flux.shape

        index  = np.arange(np.max(shape))
        pixcrd = np.column_stack([index]*len(shape))

        # 使用WCS对象将像素坐标转换为世界坐标，这里的世界坐标即为波长
        world = wcs.wcs_pix2world(pixcrd, 0)
        
        # 假设波长信息是在世界坐标的最后一个维度
        wavelength = world[:, 0]

        time=self.read_time(header,time_key,time_format)# time_key=None,time_format=None
            
        return spec_flux, wavelength, header, time
    
    class spectrum():
        def __init__(self,spec_name=None,
                          wv       =np.array([]),
                          flux     =np.array([]),
                          flux_e   =np.array([]),
                          time =np.nan,
                          header   =fits.Header(),
                          exptime  =np.nan,telescope=None,
                          flux_unit=u.Unit("erg cm-2 s-1 AA-1"),
                          wv_unit  =u.AA,
                          peak_time=Time(0,format="jd"),
                          spec_id=np.nan,
                          redshift=np.nan):
            self.spec_name=spec_name
            if ~any(flux_e):
                flux_e=np.zeros(len(wv))
            
            I=(~np.isnan(wv))*(~np.isnan(flux))*(~np.isnan(flux))
            self.wv=wv[I]
            self.flux=flux[I]
            self.flux_e=flux_e[I]
            self.time=time
            self.header=header
            self.exptime=exptime
            self.telescope=telescope
            self.peak_time=peak_time
            self.spec_id=spec_id

            self.flux_unit=flux_unit
            self.wv_unit  =wv_unit
            self.flux_u   =self.flux*flux_unit
            self.flux_e_u =self.flux_e*flux_unit
            self.wv_u     =self.wv*wv_unit
            
            self.flux_smooth_param=[]
            self.flux_smooth=None
            self.SN_name=""
            self.redshift=redshift
        @property
        def wv_flux(self):
            return np.stack((self.wv,self.flux), axis=1)
        @property
        def wv_z(self):
            if np.isnan(self.redshift):
                Warning.warn("redshift is not set.")
            elif self.redshift==0.0:
                Warning.warn("redshift is 0.0")
            return self.wv/(1+self.redshift)
        @property
        def flux_n(self):
            return self.flux/np.nanmedian(self.flux)
        def plot(self,color="k",**arg):
            # plt.errorbar()
            # spec_mjd=float(data_json[name]['spectra'][spec_id]["time"])
            print(self.peak_time==Time(0,format="jd"))
            peak_mjd = 0 if self.peak_time==Time(0,format="jd") else self.peak_time.mjd
            plt.plot(self.wv,
                     self.flux,
                     label="%s,phase=%.2f(%d)"%(self.time.isot,self.time.mjd-peak_mjd,self.spec_id),color=color,**arg)
            plt.xlabel('Wavelength[A]')
            plt.ylabel('Normalized Flux')
            plt.legend()



        def spec_save(self,):
            spec_name= []
            name     = self.SN_name
            spec_name.append(name)
            spec_name.append(self.spec_id)
            if self.time.mjd >30000: spec_name.append(self.time.to_value("isot").split("T")[0].replace("-",""))
            if "telescope" in self.__dict__: spec_name.append(self.telescope.replace(" ","").replace("-","_"))

            new_file_name="_".join(spec_name)+".dat"
            
            np.savetxt(fname=new_file_name,
                       X=np.stack((self.wv,self.wv_flux),axis=1))
            
            print(os.path.abspath(new_file_name), "is saved.")
    def read_time(self,header,key=None,fmt=None):
        from astropy.time import Time
        possible_times_fmt = {'DATE-OBS':"isot", 
                                'MJD-OBS':"mjd", 
                                'MJD':"mjd", 
                                'JD':"jd", 
                                "ISOT":"isot", 
                                "JD":"jd",}
        if any(fmt) and any(key):
            return Time(header[key],format=fmt)
        else:
            for k in possible_times_fmt.keys():
                if k in header:
                    print(k,":",header[k])
                    return Time(header[k],format=possible_times_fmt[k])
        return Time(0,format="jd")


    @property    
    def all_spec_info(self):
        import pandas as pd
        data=[[self.SpecAll[i].spec_name,i,self.SpecAll[i].time.mjd] for i in range(len(self.SpecAll))]
        return pd.DataFrame(data,columns=["name","i","mjd"])
            
    def sort_SpecTime(self,):
        spec_fits_sort=self.all_spec_info.sort_values(by="mjd").reset_index(drop=True)
        SpecAll_sort=[]
        for i in spec_fits_sort.index:
            self.SpecAll[spec_fits_sort.iloc[i]["i"]].spec_id=i
            self.SpecAll[spec_fits_sort.iloc[i]["i"]].SN_name=self.name
            SpecAll_sort.append(self.SpecAll[spec_fits_sort.iloc[i]["i"]])
        self.SpecAll=SpecAll_sort

    def add_LJ2m4(self,spec_path):
        spec_flux, wave, header, time= self.readspec(spec_path,time_key="MJD-OBS",time_format="mjd")
        new_spectrum=self.spectrum(
                                   flux=spec_flux[0,0,:],
                                   flux_e=spec_flux[1,0,:],
                                   wv=wave,
                                   spec_name=os.path.basename(spec_path),
                                   telescope="LJT",
                                   header=header,
                                   time=time,
                                   redshift=self.redshift,
                                   exptime=header["EXPTIME"])
        self.SpecAll.append(new_spectrum)
        self.sort_SpecTime()

    
    def add_XL2m16(self,spec_path):
        spec_flux, wave, header, time= self.readspec(spec_path,time_key="JD",time_format="jd")
        new_spectrum=self.spectrum(flux=spec_flux,
                                   wv=wave,
                                   spec_name=os.path.basename(spec_path),
                                   telescope="XLT",
                                   header=header,
                                   redshift=self.redshift,
                                   time=time,
                                   exptime=header["EXPTIME"])
        self.SpecAll.append(new_spectrum)
        self.sort_SpecTime()
    # Xinglong 2.16m
    def add_Lick_Kast(self,spec_path):
        spec_flux, wave, header, time= self.readspec(spec_path,time_key="DATE-OBS",time_format="isot")
        new_spectrum=self.spectrum(flux=spec_flux[0,:]*1e-16,
                                   flux_e=spec_flux[1,:]*1e-16,
                                   wv=wave,
                                   spec_name=os.path.basename(spec_path),
                                   telescope="Shane",
                                   header=header,
                                   time=time,
                                   redshift=self.redshift,
                                   exptime=header["EXPTIME"])
        self.SpecAll.append(new_spectrum)
        self.sort_SpecTime()

    def add_Lick_Kast_coadd(self,spec_name):
        spec_flm=open(spec_name).readlines()
        host_wv=[]
        host_flx=[]
        host_flxe=[]
        for i in spec_flm:
            i0=i.strip("\n").split()
            host_wv  .append(float(i0[0]))
            host_flx .append(float(i0[1]))
            host_flxe.append(float(i0[2]))
        wv  =np.array(host_wv)
        flx =np.array(host_flx)
        flxe=np.array(host_flxe)
        mjd0=spec_name.strip(".coadd.blotch-ui.flm").split("-")[-1]
        mjd=float(mjd0[-4:])+Time({"year":int(mjd0[0:4]),"month":int(mjd0[4:6]),"day":int(mjd0[6:8])}).to_value("mjd")
        header={}

        new_spectrum=self.spectrum(flux=flx,
                                   flux_e=flxe,
                                   wv=wv,
                                   spec_name=os.path.basename(spec_name),
                                   telescope="Shane",
                                   header=header,
                                   time=Time(mjd,format="mjd"),
                                   redshift=self.redshift,
                                   exptime=header["EXPTIME"])
        self.SpecAll.append(new_spectrum)
        self.sort_SpecTime()
        
    def add_Keck(self,spec_path):
        spec_txt  = open(spec_path).readlines()
        data_line = spec_txt.index(list(filter(lambda x: x[0:2]== "##",spec_txt))[0])+1
        header=fits.Header()
        for  i in spec_txt[:data_line-1]:
            if      "# COMMENT" in i[:9]:
                pass
            elif "# HISTORY" in i[:9] or "# END" in i[:5]:
                pass
            else: 
                tmp=i.strip("# ").strip("\n")
                item  =tmp.split("=")[0]
                coment=tmp.split("/ ")[-1]
                value=tmp.strip(item+"= ").strip("/ "+coment)
                item0 = item.strip(" ")
                if "\'" in value:
                    value=value.strip(value.split("\'")[-1]).strip("\'")
                header[item0]=value
                # if specname=1
                
        mjd  = float(header["MJD-OBS"])
        wv   = []
        flux = []
        fluxe= []
        for i in spec_txt[data_line:]:
            wv  .append(float(i.split()[0]))
            flux.append(float(i.split()[1]))
            fluxe.append(float(i.split()[2]))
        wv   = np.array(wv)
        flux = np.array(flux)
        header["MJD-OBS"]=mjd

        new_spectrum=self.spectrum(flux=flux,
                                   flux_e=fluxe,
                                   wv=wv,
                                   spec_name=os.path.basename(spec_path),
                                   telescope="Keck",
                                   header=header,
                                   time=Time(mjd,format="mjd"),
                                   redshift=self.redshift,
                                   exptime=header["EXPTIME"])
        self.SpecAll.append(new_spectrum)
        self.sort_SpecTime()
    
    def load_fits_spec(self,File_Name,
                       instrument,
                        time_key,time_format,
                        exptime_key):
        spec_flux, wave, header, time= self.readspec(File_Name,time_key=time_key,time_format=time_format)

        if len(spec_flux.shape)>1:
            spec_flux=spec_flux[:,0]
            spec_fluxe=spec_flux[:,1]
        else:        
            spec_fluxe=np.zeros(len(spec_flux))

        new_spectrum=self.spectrum(flux=spec_flux,
                                   flux_e=spec_fluxe,
                                   wv=wave,
                                   spec_name=os.path.basename(File_Name),
                                   telescope=instrument,
                                   header=header,
                                   time=time,
                                   redshift=self.redshift,
                                   exptime=header[exptime_key])
        self.SpecAll.append(new_spectrum)
        self.sort_SpecTime()    
    def load_txt_spec(self,File_Name,mjd,instrument="Unkown",scale=1,header=fits.Header()):
        import numpy as np
        data = np.loadtxt(File_Name)    
        new_spectrum=self.spectrum(flux=data[:,1]*scale,
                                   wv=data[:,0],
                                   spec_name=os.path.basename(File_Name),
                                   telescope=instrument,
                                   header=header,
                                   time=Time(mjd,format="mjd"),
                                   redshift=self.redshift,
                                   exptime=np.nan)
        self.SpecAll.append(new_spectrum)
        self.sort_SpecTime()
    def add_LCO(self,spec_path):
        spec_flux, wave, header, time= self.readspec(spec_path,time_key="MJD-OBS",time_format="mjd")
        new_spectrum=self.spectrum(flux=spec_flux[0,0,:],
                                   flux_e=spec_flux[1,0,:],
                                   wv=wave,
                                   spec_name=os.path.basename(spec_path),
                                   telescope=header["ORIGNAME"][:5],
                                   header=header,
                                   time=time,
                                   redshift=self.redshift,
                                   exptime=header["EXPTIME"])
        self.SpecAll.append(new_spectrum)
        self.sort_SpecTime()

    def SNeSpace(self,
                 json_path="/home/ysy/data/SNe.Space/all_json/"):
        import json
        import json,os
        import numpy as np
        from astropy.time import Time
        SN_name=self.name
        


        if os.path.exists(json_path):
            with open(json_path+SN_name+".json") as js:
                data_json= json.load(js)

            # all_band=np.array([i["band"] for i in data_json[SN_name]["photometry"] if "band" in i])
            # self.all_band={f:len(np.where(all_band==f)[0]) for f in np.unique(all_band)}
            self.data_json=data_json
            self.name=SN_name
            self.spec_num=len(data_json[SN_name]['spectra'])
            self.redshift=float(data_json[SN_name]['redshift'][0]["value"])
            
            spec_mjd_all = [float(data_json[SN_name]['spectra'][i]["time"]) for i in range(self.spec_num)]
            if any(data_json[SN_name]['maxdate']):
                peak_date=Time(data_json[SN_name]['maxdate'][0]["value"].replace("/","-"))
            else:
                peak_date=Time(np.nanmin(spec_mjd_all),format="mjd")
            
            self.ref_mjd  =peak_date.to_value("mjd")
            
            all_sources={}
            print(data_json[SN_name]["sources"])
            for sc in data_json[SN_name]["sources"]:
                if "alias" in  sc and "bibcode" in sc:
                    all_sources[sc["alias"]]=sc["bibcode"]
                elif "alias" in  sc and 'url' in sc:
                    all_sources[sc["alias"]]=sc["url"]
                else:
                    continue
            if "source" in data_json[SN_name]['redshift'][0]:
                self.redshift_bibcode=[all_sources[i]for i in data_json[SN_name]['redshift'][0]["source"].split(",")]
            time_format={'MJD':"mjd"}
            # print(data_json[SN_name]['sources'])
            # print(data_json[SN_name]['spectra'][0]['source'])
            # print(data_json[SN_name]['spectra'][0].keys())
            for i in range(self.spec_num):
                time=np.nan
                spec=data_json[SN_name]['spectra'][i]

                if "time" in spec and 'u_time' in spec:
                    time=Time(spec["time"],format=time_format[spec['u_time']])
                else:
                    continue
                
                if "source" in spec:
                    source=[all_sources[i]   for i in spec["source"].split(",") if i in all_sources]
                else:
                    continue
                
                telescope=spec["instrument"] if "instrument" in spec else "unknown"

                new_spectrum= self.spectrum(flux=np.float32(spec["data"])[:,1],
                                            wv=np.float32(spec["data"])[:,0],
                                            spec_name=os.path.basename(f"SNeSpace_{i}"),
                                            telescope=telescope,
                                            header=fits.Header(),
                                            redshift=self.redshift,
                                            time=time,
                                            )
                new_spectrum.bibcode=source
                self.SpecAll.append(new_spectrum)
                self.sort_SpecTime()        


class LightCurve():
    def __init__(self,band="",name="",redshift=0.0,distmod=0,ebv_host=[np.nan,np.nan],ebv_MW=[np.nan,np.nan],csv_name="",exptime=0,obs_mag=False):
        import pandas as pd
        self.name    =name
        self.redshift=redshift
        self.allband=[i for i in band ]
        self._distmod = distmod
        self._ebv_host= ebv_host
        self._ebv_MW  = ebv_MW

        self.A_lambda_host    ={}
        self.A_lambda_host_err={}
        self.A_lambda_MW      ={}
        self.A_lambda_MW_err  ={}
        self.A_lambda    ={}
        self.A_lambda_err={}
        self._lc_obs_pd=pd.DataFrame(columns="mjd mag mag_e detection instrument filter".split())
        self.exptime=exptime
        self._obs_mag=obs_mag
        for f in self.allband:            
            self.determine_extinction(f)
            setattr(self,f,
                    self.LightCurve_data(
                                         self,
                                         Filter      =f,
                                         A_lambda    =self.A_lambda[f],
                                         A_lambda_err=self.A_lambda_err[f],
                                            ))
            if f not in self.allband:self.allband.append(f)
        if os.path.exists(csv_name):
            # self.allband=np.unique(pd.read_csv(csv_name)["filter"])
            self.load_csv_data(csv_name)
    @property
    def json_drop(self):
        drop_list={f:{i:[] for i in np.unique(self.__dict__[f].instrument)} for f in self.allband}
        return drop_list
    def drop_list_save(self):
        import json
        with open("./%s_drop_list.json"%self.name, 'w') as json_file:
            json.dump(self.json_drop,json_file,ensure_ascii=False, indent=4)
        print("%s_drop_list.json"%self.name,"has been saved in",os.path.abspath("./"))
    def drop_data(self,json_path="",Print=True):
        import json
        if os.path.exists(json_path):
            with open(json_path, 'r') as json_file:
                self.drop_list = json.load(json_file)
        for f in self.drop_list:
            for i in self.drop_list[f]:
                if any(self.drop_list[f][i]):
                    self.__dict__[f].dropdata(self.drop_list[f][i],i,Print=Print)
    @property
    def lc_obs_pd(self):
        have_data=False
        for f in self.allband:
            if any(self.__dict__[f].mjd):
                have_data=True
                break
        if not have_data:
            return pd.DataFrame(columns="mjd mag mag_e detection instrument filter".split())
        else:
            for f in self.allband:
                self._lc_obs_pd=self._lc_obs_pd.append(self.__dict__[f].lc_obs_pd[self._lc_obs_pd.columns])
            self._lc_obs_pd=self._lc_obs_pd.reset_index(drop=True)
        return self._lc_obs_pd.drop_duplicates()
    def __getitem__(self, key):
        if key in self.__dict__.keys():
            return self.__dict__[key]
        else:
            raise KeyError(f"Key '{key}' not found")
    @property
    def lc_pd(self):
        # import pd
        All_data=[]
        for f in self.allband:
            All_data=[list(i)+[f] for i in zip( self.__dict__[f].mjd,
                                            self.__dict__[f].mag,
                                            self.__dict__[f].mag_e,
                                            self.__dict__[f].detection,
                                            self.__dict__[f].instrument)
                                            ]+All_data
        data=pd.DataFrame(data=All_data,columns="mjd mag mag_e detection instrument filter".split())
        data= data.sort_values(by='mjd')
        return data.drop_duplicates()
    
    def save_lc_pd(self,path=""):
        if path=="":
            path=os.path.abspath("./")+self.name+".csv"
        self.lc_pd.to_csv(path, index=False, mode='w')
        print("save extinction-corrected lightcurves to:",os.path.abspath(path))
    def save_lc_obs_pd(self,path=""):
        if path=="":
            path=os.path.abspath("./")+self.name+".csv"
        self.lc_obs_pd.to_csv(path, index=False, mode='w')
        print("save extinction-corrected lightcurves to:",os.path.abspath(path))
        # import pd
        
    # ["mjd","mag","mag_e","instrument","filter","detection"]

    @property
    def all_instrument(self):
        return np.unique(self.lc_pd["instrument"].values)
    
    def load_csv_data(self,csv_name):
        import pandas as pd
        lc_pd = pd.read_csv(csv_name)
        lc_pd = lc_pd[['mjd', 'mag', 'mag_e',"instrument","filter","detection"]]
        add_filters=np.unique(lc_pd["filter"].values)


        for f in add_filters:            
            self.determine_extinction(f)
            if f not in self.allband:
                setattr(self,f,
                        self.LightCurve_data(
                                            self,
                                            Filter      =f,
                                            A_lambda    =self.A_lambda[f],
                                            A_lambda_err=self.A_lambda_err[f],
                                            ))
                self.allband.append(f)

            self.__dict__[f].add_LightCurve_data(lc_pd.values.tolist())

    def determine_extinction(self,band,extco_host=3.1):
        
        if type(self._ebv_host)==float:
            self._ebv_host=[self._ebv_host,self._ebv_host*0.1]
            warnings.warn("The Error of ebv of the host galayx is set as 0.2 of ebv_Host")
        if type(self._ebv_MW  )==float:
            self._ebv_MW=[self._ebv_MW,self._ebv_MW*0.1]
            warnings.warn("The Error of ebv of the Milky Way is set as 0.2 of ebv_MW")
        
        if np.isnan(self._ebv_MW[0]):
            warnings.warn("the ebv of host Miky Way to nan. This may lead to unexpected behavior.")
        if np.isnan(self._ebv_host[0]):
            warnings.warn("the ebv of host galaxy is nan. This may lead to unexpected behavior.")
        
        self.extco_all= {
            #Extinction coefficients in A_lam / E(B-V). Uses York Extinction Solver (http://www.cadc-ccda.hia-iha.nrc-cnrc.gc.ca/community/YorkExtinctionSolver/coefficients.cgi)
                'u': 4.786,  'g': 3.587, 'r': 2.471,'C': 2.471, 'i': 1.798,  'z': 1.403, 'y': 1.228, 'w':2.762, 'Y': 1.228,
                'U': 4.744,  'B': 4.016, 'V': 3.011, 'R': 2.386, 'G': 2.216, 'I': 1.684, 'J': 0.813, 'H': 0.516,
                'K': 0.337, 'S': 8.795, 'D': 9.270, 'A': 6.432,  'F': 8.054,  'N': 8.969, 'o': 2.185, 'c': 3.111,
                'W': 0.190, 'Q': 0.127, 
                'Ks': 0.337,"u'": 4.786,  "g'": 3.587, "r'": 2.471, "i'": 1.798,
                # https://academic.oup.com/mnras/article/525/2/2701/7244724?login=false
                # "UVW2": 5.6,"UVM2": 6.99,"UVW1": 4.9,
                "UVW2": 8.795, "UVM2": 9.270, "UVW1": 6.432,
                'uuu': 4.13,"ubb": 3.41,"uvv": 4.27,
                }

        if band not in self.extco_all:
            print(self.extco_all)
            raise KeyError("The %s band extinction coeffceincy not in the list"%band)
        else:
            self.A_lambda_MW    [band]  = self.extco_all[band]*self._ebv_MW[0]
            self.A_lambda_MW_err[band]  = self.extco_all[band]*self._ebv_MW[1]
        
        # if extco_host==3.1:
        #     self.A_lambda_host    [band] =self.extco_all[band]*self._ebv_host[0]
        #     self.A_lambda_host_err[band] =self.extco_all[band]*self._ebv_host[1]
        # else:
        #     self.A_lambda_host    [band]= extco_host * self._ebv_host[0] 
        #     self.A_lambda_host_err[band]= extco_host * self._ebv_host[1]
        
        self.A_lambda_host    [band] =self.extco_all[band]*self._ebv_host[0]
        self.A_lambda_host_err[band] =self.extco_all[band]*self._ebv_host[1]
        
        self.A_lambda    [band] =self.A_lambda_host[band]+self.A_lambda_MW[band]
        self.A_lambda_err[band] =(self.A_lambda_host_err[band]**2+self.A_lambda_MW_err[band]**2)**0.5

        return self.A_lambda_host, self.A_lambda_MW, self.A_lambda,

    @property
    def distmod(self):
        if self._distmod == 0:
            if self.redshift != 0:
                from astropy.coordinates import Distance
                self._distmod=Distance(z=self.redshift).distmod.value
            else:
                warnings.warn("distmod is equal to 0. The redshift=0")

        return self._distmod
    
    @property
    def ebv_MW(self):
        if np.isnan(self._ebv_MW[0]):
            warnings.warn("the ebv of MW is equal to nan. This may lead to unexpected behavior.")
        return self._ebv_MW
    @property
    def ebv_host(self):
        if np.isnan(self._ebv_host[0]):
            warnings.warn("the ebv of host galaxy to nan. This may lead to unexpected behavior.")
        return self._ebv_host
    
    def add_band(self,band):
        if type(band)==str:
            band=list(band)
        for b in band:
            if b not in self.allband:
                self.allband.append(b)
                self.determine_extinction(b)
                setattr(self,b,
                        self.LightCurve_data(
                                             self,
                                             Filter      =b,
                                             A_lambda    =self.A_lambda[b],
                                             A_lambda_err=self.A_lambda_err[b],
                                             ))
            else:
                print(Warning("%s has already in the Light Curve"%b))
        # self.photo_num+=
    def combine_photometry(self,combined_band=""):
        if not any(combined_band):
            combined_band=self.allband
        for f in combined_band:
            self.__dict__[f].combine_data()
    
    # def color():
    #     mjd=1
    #     color=2
    #     return mjd, color
    @property
    def obs_mag(self):
        return self._obs_mag

    @obs_mag.setter
    def obs_mag(self, value):
        self._obs_mag = value
    class LightCurve_data:

        def __init__(self,outer_instance,Filter="Unknown",A_lambda=0.0,A_lambda_err=0.0):
            self._mjd       =np.array([])
            self._mag       =np.array([])
            self._mag_obs   =np.array([])
            self._mag_e     =np.array([])
            self._mag_obs_e =np.array([])
            self._data      ="Umpty"
            self._if_ExtCor =np.array([])
            self.filter     =Filter
            self.instrument =[]
            self.detection  =[]
            self.A_lambda    =A_lambda
            self.A_lambda_err=A_lambda_err
            self.peak_mjd    =np.nan
            self.Vacca_Leibundgut = {}
            self._outer_instance = outer_instance
            self._lc_obs_pd=pd.DataFrame(columns=["mjd","mag","mag_e","instrument","filter","detection"])
        
        @property
        def obs_mag(self):
            return self._outer_instance.obs_mag

        @obs_mag.setter
        def obs_mag(self, value):
            self._outer_instance.obs_mag = value
        
        @property
        def photo_num(self):
            return len(self._mjd)
        @property
        def lc_obs_pd(self):
            if self.photo_num==0:
                return pd.DataFrame(columns=["mjd","mag","mag_e","instrument","filter","detection"])
            return self._lc_obs_pd
        # def combine_data(self):
        #     self
        
        @property
        def mag(self):
            if self.obs_mag:
                return self._mag_obs
            if any(self._if_ExtCor==False):
                self.extinction_correction()
                return self._mag
            else:
                return self._mag

        @property
        def mag_e(self):
            if self.obs_mag:
                return self._mag_obs_e
            if any(self._if_ExtCor==False):
                self.extinction_correction()
                return self._mag_e
            else:
                return self._mag_e
        @property
        def mjd(self):
            return self._mjd
            
        def extinction_correction(self):
            false_indices = np.where(self._if_ExtCor == False)
            # print(self.A_lambda_err)
            self._mag  [false_indices] -= self.A_lambda
            self._mag_e[false_indices] = (self.A_lambda_err**2 + self._mag_e[false_indices]**2)**0.5
            self._if_ExtCor =np.array(len(self._mag)*[True])
                
        def filter_and_sort(self,a1, a2, a3):
            # 获取a1从小到大排序后的索引
            sorted_indices = np.argsort(a1)
            # 筛选和排序a2，a3
            valid_indices = []
            for idx in sorted_indices:
                if not np.isnan(a2[idx]) and not np.isnan(a3[idx]):
                    valid_indices.append(idx)
            return valid_indices
        
        def add_LightCurve_data(self,data_load,instrument="all"):
            LightCurve_tmp0 = list(filter(lambda X : self.filter == X[4],data_load))
            if instrument=="all":
                LightCurve_tmp1 = LightCurve_tmp0
            else:
                LightCurve_tmp1 = list(filter(lambda X : instrument in X[3],LightCurve_tmp0  ))

            import pandas as pd
        
        
            add_data=pd.DataFrame(LightCurve_tmp1,columns=["mjd","mag","mag_e","instrument","filter","detection"])
            self._lc_obs_pd=self._lc_obs_pd.append(add_data,ignore_index=True)
            self._lc_obs_pd=self._lc_obs_pd.drop_duplicates()
            self._lc_obs_pd=self._lc_obs_pd.sort_values(by=["mjd"]).reset_index(drop=True)


            for i in LightCurve_tmp1:
                self._mjd        = np.append(self._mjd        ,i[0])
                self._mag        = np.append(self._mag        ,i[1])
                self._mag_obs    = np.append(self._mag_obs    ,i[1])
                self._mag_e      = np.append(self._mag_e      ,i[2])
                self._mag_obs_e  = np.append(self._mag_obs_e  ,i[2])
                self.detection   = np.append(self.detection  ,i[5])
                self.instrument  = np.append(self.instrument ,i[3])
                self._if_ExtCor  = np.append(self._if_ExtCor ,[False])
            

            new_index=self.filter_and_sort(self._mjd,self._mag,self._mag_e)
            self._mjd        =self._mjd[new_index]
            self._mag        =self._mag[new_index]
            self._mag_e      =self._mag_e[new_index]
            self._mag_obs    =self._mag_obs[new_index]
            self._mag_obs_e  =self._mag_obs_e[new_index]
            self.detection   =self.detection[new_index]
            self.instrument  =self.instrument[new_index]
            self._if_ExtCor  =self._if_ExtCor[new_index]
            
            self.extinction_correction()
            # self.host_galaxy = np.append(self.host_galaxy,np.zeros(len(self.mjd)))
        def dropdata(self,drop_index,instrument,Print=True):
            # 初始化结果列表
            float_indices_values = []
            int_indices_values   = []
            # 遍历列表并记录索引和值
            for index, number in enumerate(drop_index):
                if isinstance(number, float) or isinstance(number, np.float64):
                    float_indices_values.append((index, number))
                elif isinstance(number, int) or isinstance(number,np.int0):
                    int_indices_values.append((index, number))
            ins=self.instrument==instrument
            all_drop_indice=[]
            if any(int_indices_values):
                for i in int_indices_values:
                    if i[1] < len(np.where(ins)[0]):
                        indice=np.where((self._mjd==self._mjd[ins][i[1]])*(self.instrument==instrument))[0][0]
                        all_drop_indice.append(indice)
            # print(float_indices_values)
            if any(float_indices_values):
                for d in float_indices_values:
                    if d[1] in self.mjd[ins]:
                        indice=np.where((self._mjd==d[1])*(self.instrument==instrument))[0][0]
                        all_drop_indice.append(indice)
            drp_i=np.array(all_drop_indice)
            if Print:
                print(instrument,self.filter)            
                print(self._mjd[drp_i])
                print(self._mag[drp_i])
                print(self. mag[drp_i])
                print(self.mag_e[drp_i])
                print(self.detection[drp_i])
                print(self._if_ExtCor[drp_i])
            self._mjd        = np.delete(self._mjd      ,drp_i)
            self._mag        = np.delete(self._mag      ,drp_i)
            self._mag_obs    = np.delete(self._mag_obs  ,drp_i)
            self._mag_e      = np.delete(self._mag_e    ,drp_i)
            self._mag_obs_e  = np.delete(self._mag_obs_e,drp_i)
            self.detection   = np.delete(self.detection ,drp_i)
            self.instrument  = np.delete(self.instrument,drp_i)
            self._if_ExtCor  = np.delete(self._if_ExtCor,drp_i)
            self._lc_obs_pd  = self._lc_obs_pd.drop(drp_i).reset_index(drop=True)
def LC2latex_2(LightCurve,Filter=[],inst="",ref_phase=0):
    import pandas as pd
    
    data_=LightCurve.lc_obs_pd["mjd mag mag_e detection instrument filter".split()]
    All_instrument=LightCurve.all_instrument

    if inst=="":
        print_inst=All_instrument
    else:
        print_inst=[ins for ins in All_instrument if ins in  inst]
    
    data_=data_[data_["instrument"].isin(print_inst)].reset_index(drop=True)
    
    if not any(Filter):    
        print_band=np.unique(data_["filter"].values)
    else:
        print_band=[f for f in np.unique(data_["filter"].values) if f in Filter]
    time_gap=0.2
    ALL_index=[]
    for i in data_.index:
        ii=list(np.where(((data_["mjd"].values>data_["mjd"].values[i]-time_gap)&
                          (data_["mjd"].values<data_["mjd"].values[i]+time_gap)))[0])
        if ii not in ALL_index: 
            ALL_index.append(ii)
    
    INSTRU_name={
        "Atlas":"ATLAS",
        "tnt80cm":"TNT",
        "PSP+BIN2":"NEXT",
        "HMT":"HMT",
        "ZTF":"ZTF",
        "LCO_ELP1m06":"LCO_ELP1m06",
        "LCO_ELP1m08":"LCO_ELP1m08",
        "LCO_tfn0m414":"LCO_tfn0m414",
        "LCO_tfn1m01":"LCO_tfn1m01",   
    }
    # ALL_print_data=[]
    new_column= ["MJD","Phase","Instrument"]
    for b in print_band:
        new_column += [b,b+"_err",b+"_det"]
                                        
    ALL_print_data=pd.DataFrame(columns=new_column)
    for sd in ALL_index:
        for inst in print_inst:
            data_inst_day = data_.iloc[sd][data_["instrument"]==inst]
            # print(data_inst_day)
            data_inst_day = data_inst_day[data_inst_day['filter'].isin(print_band)]
            if data_inst_day.values.any():
                new_row=pd.Series([np.nan] * len(ALL_print_data.columns),
                                           index=ALL_print_data.columns)
                # mag_sd_sf={}
                for f in print_band:
                    if 0.0 in  data_inst_day["detection"] and 1.0 in data_inst_day["detection"]:
                        # data_inst_day.drop((data_inst_day["detection"]==0.0).index).reset_index(drop=True)
                        data_inst_day=data_inst_day[data_inst_day["detection"]!=0.0].reset_index(drop=True)
                    
                    new_row["MJD"]    =np.mean(data_inst_day["mjd"      ].values)
                    new_row["Phase"]  =np.mean(data_inst_day["mjd"      ].values)-ref_phase
                    new_row[f]        =np.mean(data_inst_day["mag"      ][data_inst_day["filter"]==f].values)
                    new_row[f+"_err"] =np.mean(data_inst_day["mag_e"    ][data_inst_day["filter"]==f].values**2)**0.5
                    new_row[f+"_det"] =np.mean(data_inst_day["detection"][data_inst_day["filter"]==f].values)
                    new_row["Instrument"]=INSTRU_name[inst] if inst in INSTRU_name else inst
                ALL_print_data = ALL_print_data.append(new_row, 
                                                        ignore_index=True)
    ALL_print_data.sort_values(by='MJD').reset_index(drop=True)    
    

    All_item="MJD & Phase & Instrument"
    print_fmt="%.2f & %.2f & %s "
    for f in print_band:
        All_item += " & "
        All_item += f
    print("\\hline")
    print(All_item+"\\\\")
    print("\\hline")
    print("\\hline")
    
    for i in ALL_print_data.index:
        MAG_str=""
        for f in print_band:
            # 无数据
            if np.isnan(ALL_print_data.iloc[i][f]):
                mag_str="..."
            # 上限
            elif ALL_print_data.iloc[i][f+"_det"]==0:
                mag_str="$>$%.2f"%ALL_print_data.iloc[i][f]
            # 数据
            else:
                mag_str="%.2f$\pm$%.2f"%(ALL_print_data.iloc[i][f],
                                         ALL_print_data.iloc[i][f+"_err"])
            MAG_str+=" & "
            MAG_str+=mag_str
        # print(ALL_print_data[i])
        print(print_fmt%(
                        ALL_print_data.iloc[i]["MJD"],
                        ALL_print_data.iloc[i]["Phase"],
                        ALL_print_data.iloc[i]["Instrument"])+
                        MAG_str+
                        " \\\\")
    print("\\hline")

def LC2latex_1(LightCurve,Filter=[],inst="",ref_phase=0):
# def LC2latex(LightCurve,Filter=[],ref_phase=0):
    print(r"{ccccc}")
    print("\\hline")
    print("\\hline")
    print("MJD & Phase & Filter & Magnitude &Instrument\\\\")
    print("\\hline")

    
    INSTRU_name={
        "Atlas":"ATLAS",
        "tnt80cm":"TNT",
        "PSP+BIN2":"NEXT",
        "HMT":"HMT",
        "ZTF":"ZTF",
        "LCO_ELP1m06":"LCO_ELP1m06",
        "LCO_ELP1m08":"LCO_ELP1m08",
        "LCO_tfn0m414":"LCO_tfn0m414",
        "LCO_tfn1m01":"LCO_tfn1m01",   
    }
    All_data=LightCurve.lc_pd
    All_instrument=LightCurve.all_instrument

    if inst=="":
        print_inst=All_instrument
    else:
        print_inst=[ins for ins in All_instrument if ins in  inst]
    
    data_=All_data[All_data["instrument"].isin(print_inst)].reset_index(drop=True)
    
    if not any(Filter):
        print_band=np.unique(data_["filter"].values)
    else:
        print_band=[f for f in np.unique(data_["filter"].values) if f in Filter]
    
    #mjd, mag, mage, detection, instrument_name, filter

    for ins in print_inst:
    # for ins in INSTRU_name:
        data_inst=data_[data_["instrument"]==ins]
        
        for f in print_band:
            data_inst_filter=data_inst[data_inst["filter"]==f].sort_values(by="mjd").reset_index(drop=True)
            
            if not data_inst_filter.empty:
                
            
                phase_str=""
                for i in data_inst_filter.index:
                    data=data_inst_filter.iloc[i]

                    phase_str="%.2f"%(data["mjd"]-ref_phase)
                    if data["mjd"]-ref_phase>=0:phase_str="+"+phase_str
                    
                    
                    if data["detection"]==1:
                        print("%.2f & %s & $%s$ & %.2f$\pm$%.2f & %s \\\\"%(
                            data["mjd"],
                            phase_str,
                            f,
                            data["mag"],
                            data["mag_e"],
                            INSTRU_name[ins]))
                    else:
                        print("%.2f & %s & $%s$ & $>$%.2f & %s \\\\"%(
                            data["mjd"],
                            phase_str,
                            f,
                            data["mag"],
                            INSTRU_name[ins]))
    print("\\hline")
def superbol_lc2txt(LightCurve,ref_date,path,Filter=[],inst="",mjd=[0,np.inf]):
    import pandas as pd
    superbol_filter="u g r i z B V U UVW1 UVW2 UVM2 R I".split()

    data_=LightCurve.lc_pd["mjd mag mag_e detection instrument filter".split()]
    All_instrument=LightCurve.all_instrument

    # print(All_instrument)
    # print(any((data_["detection"]==0.0).index))
    data_=data_[data_["detection"]!=0.0].reset_index(drop=True)
    data_=data_[(mjd[0]<data_["mjd"])&(data_["mjd"]<mjd[1])].reset_index(drop=True)
    if inst=="":
        print_inst=All_instrument
    else:
        print_inst=[ins for ins in All_instrument if ins in  inst]
    
    data_=data_[data_["instrument"].isin(print_inst)].reset_index(drop=True)
    
    
    if not any(Filter):
        print_band=[f for f in LightCurve.allband if f in superbol_filter]
    else:
        print_band=[f for f in LightCurve.allband if (f in Filter and f in superbol_filter)]

    time_gap=0.2
    ALL_index=[]
    for i in data_.index:
        ii=list(np.where(((data_["mjd"].values>data_["mjd"].values[i]-time_gap)&
                          (data_["mjd"].values<data_["mjd"].values[i]+time_gap)))[0])
        if ii not in ALL_index: 
            ALL_index.append(ii)
    
    INSTRU_name={
        "Atlas":"ATLAS",
        "tnt80cm":"TNT",
        "PSP+BIN2":"NEXT",
        "HMT":"HMT",
        "ZTF":"ZTF",
        "LCO_ELP1m06":"LCO_ELP1m06",
        "LCO_ELP1m08":"LCO_ELP1m08",
        "LCO_tfn0m414":"LCO_tfn0m414",
        "LCO_tfn1m01":"LCO_tfn1m01",   
    }
    # ALL_print_data=[]
    new_column= ["MJD","Instrument"]
    for b in print_band:
        new_column += [b,b+"_err"]
                                        
    ALL_print_data=pd.DataFrame(columns=new_column)
    for sd in ALL_index:
        for inst in print_inst:
            data_inst_day = data_.iloc[sd][data_["instrument"]==inst]
            # print(data_inst_day)
            data_inst_day = data_inst_day[data_inst_day['filter'].isin(print_band)]
            if data_inst_day.values.any():
                new_row=pd.Series([np.nan] * len(ALL_print_data.columns),
                                           index=ALL_print_data.columns)
                # mag_sd_sf={}
                for f in print_band:

                    new_row["MJD"]    =np.mean(data_inst_day["mjd"      ].values)
                    new_row[f]        =np.mean(data_inst_day["mag"      ][data_inst_day["filter"]==f].values)
                    new_row[f+"_err"] =np.mean(data_inst_day["mag_e"    ][data_inst_day["filter"]==f].values**2)**0.5
                    new_row["Instrument"]=inst
                ALL_print_data = ALL_print_data.append(new_row, 
                                                        ignore_index=True)
                
    header_fit_band=[]
    for i in print_band: header_fit_band=header_fit_band+[i,"\t"]
    header_filter="\t".join(["#MJD"]+header_fit_band)
    import sys
    original_stdout = sys.stdout

    with open(path,"w") as f:
        sys.stdout = f
        print(header_filter)
        for inst in print_inst:
            print("\n")
            print("#%s"%inst,)
            data_inst=ALL_print_data[ALL_print_data["Instrument"]==inst].reset_index(drop=True)
            data_inst=data_inst.sort_values(by="MJD").reset_index(drop=True)
            # print(data_inst)
            for d in data_inst.index:
                MAG_str="%.2f"%(data_inst.iloc[d]["MJD"]-ref_date)
                for f in print_band:
                    # 无数据
                    if np.isnan(data_inst.iloc[d][f]):
                        mag_str="n\tn"
                    # 数据
                    else:
                        mag_str="%.2f\t%.2f"%(data_inst.iloc[d][f],
                                            data_inst.iloc[d][f+"_err"])
                    MAG_str+="\t"
                    MAG_str+=mag_str
                print(MAG_str)
        sys.stdout = sys.__stdout__ 
    sys.stdout = original_stdout
    return ALL_print_data
def read_ZTF(path_det="",path_non=""):
    import pandas as pd
    ZTF_detection    =pd.read_csv(path_det)
    ZTF_non_detection=pd.read_csv(path_non)

    ZTF_r=[list(i) for i in zip(
        ZTF_detection[ZTF_detection['fid']==2]["mjd"].array,
        ZTF_detection[ZTF_detection['fid']==2]["magpsf"].array,
        ZTF_detection[ZTF_detection['fid']==2]["sigmagapbig"].array,
    len(ZTF_detection[ZTF_detection['fid']==2]["mjd"])*["ZTF"],
    len(ZTF_detection[ZTF_detection['fid']==2]["mjd"])*["r"],
    len(ZTF_detection[ZTF_detection['fid']==2]["mjd"])*[1],
    )]

    ZTF_g=[list(i) for i in zip(
        ZTF_detection[ZTF_detection['fid']==1]["mjd"].array,
        ZTF_detection[ZTF_detection['fid']==1]["magpsf"].array,
        ZTF_detection[ZTF_detection['fid']==1]["sigmagapbig"].array,
    len(ZTF_detection[ZTF_detection['fid']==1]["mjd"])*["ZTF"],
    len(ZTF_detection[ZTF_detection['fid']==1]["mjd"])*["g"],
    len(ZTF_detection[ZTF_detection['fid']==1]["mjd"])*[1],
    )]
    if not any(path_non):
        return ZTF_r,ZTF_g
    ZTF_r_non=[list(i) for i in zip(
        ZTF_non_detection[ZTF_non_detection['fid']==1]["mjd"].array,
        ZTF_non_detection[ZTF_non_detection['fid']==1]["diffmaglim"].array,
    len(ZTF_non_detection[ZTF_non_detection['fid']==1]["mjd"])*[1],
    len(ZTF_non_detection[ZTF_non_detection['fid']==1]["mjd"])*["ZTF"],
    len(ZTF_non_detection[ZTF_non_detection['fid']==1]["mjd"])*["r"],
    len(ZTF_non_detection[ZTF_non_detection['fid']==1]["mjd"])*[0],
    )]

    ZTF_g_non=[list(i) for i in zip(
        ZTF_non_detection[ZTF_non_detection['fid']==2]["mjd"].array,
        ZTF_non_detection[ZTF_non_detection['fid']==2]["diffmaglim"].array,
    len(ZTF_non_detection[ZTF_non_detection['fid']==2]["mjd"])*[1],
    len(ZTF_non_detection[ZTF_non_detection['fid']==2]["mjd"])*["ZTF"],
    len(ZTF_non_detection[ZTF_non_detection['fid']==2]["mjd"])*["g"],
    len(ZTF_non_detection[ZTF_non_detection['fid']==2]["mjd"])*[0],
    )]

    
    return ZTF_g+ZTF_g_non+ZTF_r+ZTF_r_non

spec_color_telescp = {
    "XL2m16": "b",
    "LJ2m4" : "g",
    "ogg2m" : "r",
    "Keck"  : "cyan",
    "Lick"  : "orange",
    "Ekar_AFOSC":"r",
    "Kast"  : "orange",
    "Unkown":"k",
    "NOT"   :"k",
    "LT"   :"coral",
    "HCT"  :"purple",
    "Other":"magenta",
    "P60"  :"olive", 
}

def extinction_MW(extinction_file):
    import pandas as pd    
    extinction_csv=pd.read_csv(extinction_file)
    A_B=extinction_csv[extinction_csv['Bandpass'].str.contains("Landolt B")]
    A_V=extinction_csv[extinction_csv['Bandpass'].str.contains("Landolt V")]
    A_B=A_B[A_B['Refcode of the publications'].str.contains("2011ApJ...737..103S")]["The Galactic extinction"].values[0]
    A_V=A_V[A_V['Refcode of the publications'].str.contains("2011ApJ...737..103S")]["The Galactic extinction"].values[0]
    
    return float(A_B-A_V)
class PlotLC():
    offset_by_filter={
    "z"  : -0.5,
    "U"  : 2.5,
    "B"  : 2.0,
    "g"  : 1.5,
    "V"  : 1.0,
    "r"  : 0.5,
    "o"  : -0.5,
    "i"  : 0.0,   
    "R"  : 0.5}
    # filled_markers = ('o', 'v', '^', '<', '>', '8', 's', 'p', '*', 'h', 'H', 'D', 'd', 'P', 'X')
    marker_by_telecope={
    "LJ2m4" :"o" ,
    "LJT" :"o" ,
    "XL2m16":"^",
    "tnt80cm":"8",
    "ELP1m" :"h",
    'LCO_ELP1m':"h",
    "TFN4dm":"*",
    "ZTF"   :"s",
    "PSP+BIN2":"p",
    "Atlas" :"X",
    "Unknown":"v",
    "N89":"D",
    "HMT":"d",
    "GOTO":"D",
    "T31":"p",
    "YaoAn":"s",
    "Swift":"*",
    }

    color_by_filter={
    'u': 'dodgerblue', 
    'g': 'g', 
    'r': 'r', 
    'i': 'magenta', 
    'z': 'k', 
    'y': '0.5', 
    'w': 'firebrick',    
    'Y': '0.5', 
    'U': 'slateblue', 
    'B': 'b', 
    'V': 'yellowgreen', 
    'R': 'crimson', 
    'G': 'salmon',
    'I': 'chocolate', 
    'J': 'darkred', 
    'H': 'orangered', 
    'K': 'saddlebrown',
    'Ks': 'saddlebrown',
    'S': 'mediumorchid', 
    'A': 'midnightblue',
    'F': 'hotpink', 
    'N': 'coral', 
    'o': 'darkorange', 
    'c': 'cyan',
    'W': 'forestgreen', 
    'Q': 'peru',
    "UVW1":"purple",
    "UVW2":"deeppink",
    "UVM2":"deepskyblue",
    "C": 'r',
    "u'": 'dodgerblue', 
    "g'": 'g', 
    "r'": 'r', 
    "z'": 'k', 
    "i'": 'magenta', 
    }
    def __init__(self,LightCurve,Xshift=0,color="default",ms=3,SizeUpLimit=(2,0.2)):
        self.LightCurve = LightCurve
        self.allband    = LightCurve.allband
        self.color      = "default"
        self.Xshift     = Xshift
        self.color      = color
        self.ms         = ms
        self.inst_plot=[]
        self.legend_handles=[]
        self.offset_filter={}
        self.plot_marker={}
        self.plot_color={}
        # self.instrument_legend=False
        
        all_marker=[i for i in plt.Line2D.markers]
        for m in [None,"None",",",0,1,2,3,4,7,8,9,11," ",""]:
            all_marker.remove(m)

        self.all_marker=all_marker
        
        self.SizeUpLimit=SizeUpLimit
        all_inst=[]
        for b in self.allband:
            for i in np.unique(self.LightCurve.__dict__[b].instrument):
                if i not in all_inst:all_inst.append(i)
        self.all_inst=all_inst
        self.marker_by_telecope={t:self.marker_by_telecope[t] for t in self.marker_by_telecope.keys()}
        # print("The fliters: ",self.allband)

                
    def drawUpLimit(self,x, y, clr,size):
        size_x=size[0]
        size_y=size[1]
        alpha=0.5
        plt.plot([x-size_x/2,x+size_x/2],[y]*2,color=clr,alpha=alpha)
        plt.plot([x]*2,[y,y+size_y],           color=clr,alpha=alpha)
        plt.scatter(x,y+size_y  ,marker="v",   color=clr,alpha=alpha)
    
    def plot_band_inst(self,band,color="default",Yshift=0,labels="default",marker="",ms=3,alpha=1,offset=0.5,
                       linestyle="None",bbox_to_anchor=(0.5, 1.15,),ncol_legend=3,draw_uplm=True):
        self.alpha      = alpha
        if any(band):
            AllPlotBand=list(band)
        else:
            AllPlotBand=self.allband

        offset_plot=Yshift
        plot_line=False
        if "-" in linestyle:
            linestyle=linestyle.replace("-","")
            plot_line=True

        for b in AllPlotBand:
            # 区分探测极限和实际探测点
            detect_index     =np.where(self.LightCurve.__dict__[b].detection==1)
            non_detect_index =np.where(self.LightCurve.__dict__[b].detection!=1)

            f      = self.LightCurve.__dict__[b].filter
            ins    = self.LightCurve.__dict__[b].instrument
            
            mjd_1  = self.LightCurve.__dict__[b].mjd  [detect_index]
            mag_1  = self.LightCurve.__dict__[b].mag  [detect_index]
            mage_1 = self.LightCurve.__dict__[b].mag_e[detect_index]

            mjd_0  = self.LightCurve.__dict__[b].mjd  [non_detect_index]
            mag_0  = self.LightCurve.__dict__[b].mag  [non_detect_index]
            mage_0 = self.LightCurve.__dict__[b].mag_e[non_detect_index]

            
            # 分配颜色
            if self.color == "default":
                if f in self.color_by_filter.keys():
                    color=self.color_by_filter[f]
                else:
                    def random_color_name():
                        import matplotlib.colors as mcolors
                        import random
                        color_names = list(mcolors.CSS4_COLORS.keys())
                        return random.choice(color_names)
                    color=random_color_name()
            if f=="C":
                facecolor="white"
            else:
                facecolor=color    

            ins_all_for_b=np.unique(ins[detect_index])
            
            # if len(ins_all)==1:ins_all=np.array(ins_all)
            label=None
            self.offset_filter[b]=offset_plot
            if plot_line:
                plt.plot(mjd_1 - self.Xshift,
                         mag_1 - offset_plot,
                         color=color, )
            self.plot_color[b]=color
            for i in ins_all_for_b:
                if any(marker):
                    pass
                elif i in self.marker_by_telecope.keys():
                    marker=self.marker_by_telecope[i]
                else:
                    m_not_in_marker_by_telecope=[m for m in self.all_marker if m not in self.marker_by_telecope.values()]
                    marker=m_not_in_marker_by_telecope[np.random.randint(0,len(m_not_in_marker_by_telecope) , 1)[0]]
                    self.marker_by_telecope[i]=marker
                
                self.plot_marker[i]=marker
                
                # if self.instrument_legend:
                    
                # else:
                #     label=labels
                self.plot_marker
                ins_index = np.where(ins[detect_index]==i)
                plt.errorbar(mjd_1[ins_index] - self.Xshift,
                            mag_1[ins_index] -offset_plot,
                            yerr  =mage_1[ins_index], 
                            color =color, 
                            ms    =ms,
                            linestyle=linestyle,
                            marker   =marker,
                            markerfacecolor=facecolor,
                            markeredgecolor=color, 
                            markeredgewidth=1, 
                            capsize=1,
                            elinewidth=1.5, 
                            capthick=2, 
                            zorder=10,
                            alpha=self.alpha,
                            label=None if i in self.inst_plot else i,
                            )
                
                # if :
                if i not in self.inst_plot: 
                    self.inst_plot.append(i)
                    self.legend_handles.append(plt.Line2D([0], [0], marker=marker, color='black', markerfacecolor="w", markersize=8, label=i))
                marker=""
                # 画上限
            # print(all_legend)
            # print(offset_plot)
            offset_plot+=offset
            if draw_uplm:
                for x,y in zip(mjd_0,mag_0):
                    self.drawUpLimit(x-self.Xshift,y+Yshift,color,size=self.SizeUpLimit)
        # if self.instrument_legend:
        legend=plt.legend(handles=self.legend_handles, handlelength=0, loc='upper center', bbox_to_anchor=bbox_to_anchor,
                                ncol=ncol_legend,columnspacing=3.0)
        legend.get_frame().set_linewidth(0)

        # print(self.marker_by_telecope)
    
        if not plt.gca().yaxis_inverted():plt.gca().invert_yaxis()
        rows, cols= plt.gca().get_subplotspec().get_gridspec().get_geometry()
        num       = plt.gca().get_subplotspec().get_topmost_subplotspec().num1
        plt.subplot(rows, cols, num+1).yaxis.set_minor_locator(plt.MultipleLocator(0.5))
        plt.subplot(rows, cols, num+1).yaxis.set_major_locator(plt.MultipleLocator(1))

        return self.plot_marker,self.plot_color

    def plot_band(self,band,color="default",Yshift=0,labels="default",marker=".",ms=3,alpha=1,linestyle="None",zorder=1,draw_uplm=True,set_ticks=True):
        self.alpha      = alpha

        if any(band):
            if type(band)==str:
                if   " " in band: AllPlotBand=band.split()
                elif "," in band: AllPlotBand=band.split(",")
                else:             AllPlotBand=list(band)
        
        else:
            AllPlotBand=self.allband
        if band in self.color_by_filter:
            AllPlotBand=[band]
        for b in AllPlotBand:
            # 区分探测极限和实际探测点
            detect_index     =np.where(self.LightCurve.__dict__[b].detection==1)
            non_detect_index =np.where(self.LightCurve.__dict__[b].detection!=1)

            f      = self.LightCurve.__dict__[b].filter
            ins    = self.LightCurve.__dict__[b].instrument
            
            mjd_1  = self.LightCurve.__dict__[b].mjd  [detect_index]
            mag_1  = self.LightCurve.__dict__[b].mag  [detect_index]
            mage_1 = self.LightCurve.__dict__[b].mag_e[detect_index]

            mjd_0  = self.LightCurve.__dict__[b].mjd  [non_detect_index]
            mag_0  = self.LightCurve.__dict__[b].mag  [non_detect_index]
            mage_0 = self.LightCurve.__dict__[b].mag_e[non_detect_index]

            # 分配颜色
            if color == "default":
                color_plot=self.color_by_filter[f]
            else:
                color_plot=color
            if f=="C":
                facecolor="white"
            else:
                facecolor=color_plot  

            ins_all_for_b=np.unique(ins[detect_index])
            all_marker=[i for i in plt.Line2D.markers]
            for m in [None,"None",",",0,1,2,3,4,7,8,9,11," ",""]:
                all_marker.remove(m)

            self.all_marker=all_marker

            # if len(ins_all)==1:ins_all=np.array(ins_all)
            
            if len(np.where(np.isnan(mag_1))[0])==len(mag_1):continue
            if np.isnan(     Yshift):continue
            if np.isnan(self.Xshift):continue
            
            # ins_index = np.where(ins[detect_index]==i)
            plt.errorbar(mjd_1 - self.Xshift,
                        mag_1 + Yshift,
                        yerr  =mage_1, 
                        color =color_plot, 
                        ms    =ms,
                        linestyle=linestyle,
                        marker   =marker,
                        markerfacecolor=facecolor,
                        markeredgecolor=color_plot, 
                        markeredgewidth=1, 
                        capsize=1,
                        elinewidth=1.5, 
                        capthick=2, 
                        zorder=zorder,
                        alpha=self.alpha,
                        label=labels,
                        )
            

            # if :
                # 画上限
            # print(all_legend)
            if draw_uplm:
                for x,y in zip(mjd_0,mag_0):
                    self.drawUpLimit(x-self.Xshift,y+Yshift,color_plot,size=self.SizeUpLimit)

    
        plt.rcParams['xtick.direction'] = 'in'
        plt.rcParams['ytick.direction'] = 'in'
        
        if not plt.gca().yaxis_inverted():plt.gca().invert_yaxis()

        ylim1,ylim2=plt.ylim()
        if set_ticks:
            if mag_1.any():
                if np.abs(ylim1-ylim2)<np.abs(np.nanmax(mag_1)-np.nanmin(mag_1))*20:
                    rows, cols= plt.gca().get_subplotspec().get_gridspec().get_geometry()
                    num       = plt.gca().get_subplotspec().get_topmost_subplotspec().num1
                    plt.subplot(rows, cols, num+1).yaxis.set_minor_locator(plt.MultipleLocator(0.5))
                    plt.subplot(rows, cols, num+1).yaxis.set_major_locator(plt.MultipleLocator(1))
    def plot_offset(self,fontsize=13,Xgap=0.01,Ygap=0.3,Xtextgap=0.01,Ytextgap=0.1,rectanglesizeX=0.04,rectanglesizeY=0.2,lable_gap=0.55):
        k=0
        self.offset_filter = dict(sorted(self.offset_filter.items(), key=lambda x: x[1],reverse=True))
        for f in self.offset_filter:
            offset_str=f
            if self.offset_filter[f]==0:
                offset_str+=""
            elif self.offset_filter[f]<0:
                offset_str+="+%.1f"%(-self.offset_filter[f])
            elif self.offset_filter[f]>0:
                offset_str+="%.1f"%(-self.offset_filter[f])
            
            yshift_plot=np.min(plt.ylim())+Ygap
            xshift_plot=np.min(plt.xlim())+Xgap
            plt.text(xshift_plot+rectanglesizeX+Xtextgap,  Ytextgap+yshift_plot+k*lable_gap, s=offset_str, fontsize=fontsize, verticalalignment='center')
            plt.gca().add_patch(plt.Rectangle((xshift_plot, yshift_plot+k*lable_gap),rectanglesizeX ,rectanglesizeY, 
                                                color=self.plot_color[f],
                                                fill=True if f !="C" else False  ,
                                                alpha=1))
            k+=1
            
def read_Swift(json_path):
    import json 
    with open(json_path,"r") as uvot_js_text:
        uvot_sub_js=json.load(uvot_js_text)
    uvot_sub=pd.DataFrame(columns=list(uvot_sub_js[0].keys()))
    for i in uvot_sub_js:
        uvot_sub=uvot_sub.append(i,ignore_index=True)
    uvot_sub["instrument"]=len(uvot_sub)*["Swift"]
    uvot_sub["detection"]=len(uvot_sub)*[0]
    uvot_sub.loc[uvot_sub['upper_limit']<uvot_sub['mag'],"detection"]=1
    uvot_sub.rename(columns={"mag_err":"mag_e"},inplace=True)
    uvot_sub = uvot_sub.drop_duplicates(subset='mjd', keep="last")
    uvot_sub=uvot_sub[["mjd","mag","mag_e","instrument","filter","detection"]]
    return uvot_sub.values.tolist()

def read_autophot_LC(csv_name):
    import numpy as np
    import pandas as pd
    REDUCED_csv=pd.read_csv(csv_name)
    num_photo_data=len(REDUCED_csv)
    all_filter=["r","i","u","g","z","B","V","I","U","C","R"]
    all_band=list(filter(lambda x: x in all_filter,REDUCED_csv.keys()))
    autophot_result=[]
    all_fits_name=[]
    
    for i in range(num_photo_data):
        reduce_csv = REDUCED_csv.iloc[i,:]
        fits_name  = reduce_csv["fname"].split("/")[-1]
        # 是否减模板
        subtract = reduce_csv.subtraction
        # 判断滤光片
        for j in all_band:
            if type(reduce_csv[j])==type(str()):
                Filter=j
                break
        mjd  = reduce_csv.mjd
        telescope=reduce_csv.telescope
        mag  = float(reduce_csv[j].strip("[").strip("]"))
        mag_e= reduce_csv[j+"_err"]
        detection_limit_SNR=3

        # if :
            # detection=0
        
        # print(reduce_csv["lmag_inject"])

        # if reduce_csv["lmag_inject"]
        
        if np.isnan(mag) and ~np.isnan(reduce_csv["lmag_inject"]):
            detection=0
            mag=reduce_csv["lmag_inject"]

        if np.isnan(reduce_csv["lmag_inject"]):
            detection= 1
        # elif reduce_csv["SNR"] > detection_limit_SNR:
        elif reduce_csv["lmag_inject"] >mag:
            detection = 1
        else:
            detection = 0
            mag       = reduce_csv["lmag_inject"]
        if  ~np.isnan(mag):
            if fits_name in all_fits_name:
                remove_item=list(filter(lambda x:np.abs(x[0]-mjd)<2e-4,autophot_result))[0]
                autophot_result.remove(remove_item)
            all_fits_name.append(fits_name)
            autophot_result.append([mjd,mag,mag_e,telescope,Filter,detection])
    return autophot_result
def read_atlas_0(txt_name):
    atlas_data=[i.split() for i in open(txt_name,"r").readlines()[1:]]
    atlas_data=[[float(z) if z!="o" else "o" for z in i[:-1]]+[i[-1]] for i in atlas_data]
    data_all=[]
    for  d in atlas_data: 
        if d[1]>0:
            detection= 1 if d[1]<d[16] else 0
            data_all=data_all+[[d[0],d[1],d[2],"ATLAS",d[5],detection]]
    return data_all




# Mag2flux=MagFlux()
# # http://svo2.cab.inta-csic.es/theory/fps/index.php?id=Swift/UVOT.UVW2_trn&&mode=browse&gname=Swift&gname2=UVOT#filter
# Mag2flux.add_filter_info({'filter':"UVW2", 
#                           'vega.flux.lambda':5.37469, 
#                           'ab.flux.lambda':25.9051, 
#                           'center_wv':2049.89,
#                           'defualt_system':"vega"})
# # http://svo2.cab.inta-csic.es/theory/fps/index.php?id=Swift/UVOT.UVW1_trn&&mode=browse&gname=Swift&gname2=UVOT#filter
# Mag2flux.add_filter_info({'filter':"UVW1", 
#                           'vega.flux.lambda':4.02204, 
#                           'ab.flux.lambda':15.7569, 
#                           'center_wv':2628.38,
#                           'defualt_system':"vega"})
# # http://svo2.cab.inta-csic.es/theory/fps/index.php?id=Swift/UVOT.UVM2_trn&&mode=browse&gname=Swift&gname2=UVOT#filter
# Mag2flux.add_filter_info({'filter':"UVM2", 
#                           'vega.flux.lambda':4.66117, 
#                           'ab.flux.lambda':21.5291, 
#                           'center_wv':2248.59,
#                           'defualt_system':"vega"})
# # print(Mag2flux.filter_info)


# # from astroquery.svo_fps import SvoFps
# # filterlist=SvoFps.get_filter_list(facility='ATLAS')
# # print(filterlist[filterlist["filterID"]=="Misc/Atlas.cyan"  ][["WavelengthRef","WidthEff"]].values)
# # print(filterlist[filterlist["filterID"]=="Misc/Atlas.orange"][["WavelengthRef","WidthEff"]].values)


# # filterlist=SvoFps.get_filter_list(facility='SLOAN')
# # print(filterlist[filterlist["filterID"]=="SLOAN/SDSS.uprime_filter"][["WavelengthRef","WidthEff"]].values)
# # filter_info=filterlist[filterlist["filterID"]=="SLOAN/SDSS.uprime_filter"].to_pandas()

# # print(filter_info.columns)
# # filter_info[[ 'filterID',"WidthEff",'ZeroPointType','ZeroPointUnit','MagSys','ZeroPoint', 'PhotSystem','TrasmissionCurve']]
# # print(filter_info.loc[0])

# from astroquery.svo_fps import SvoFps
# ilterID="SLOAN/SDSS.uprime_filter"
# ilterID="SLOAN/SDSS.zprime_filter"
# filterID='2MASS/2MASS.H'
# filterID='2MASS/2MASS.Ks'
# filterID='2MASS/2MASS.J'
# # ilterID=

# TransCurve=SvoFps.get_transmission_data(filter_id=filterID).to_pandas()
# import pyphot as pph 
# from pyphot import (unit, Filter)
# lib = pph.get_library()
# newfilter = Filter(TransCurve["Wavelength"]* unit['AA'] ,TransCurve["Transmission"], name='newfilter', dtype='photon', unit='Angstrom')

# Mag2flux.add_filter_info({'filter':"J", 
#                     'vega.flux.lambda':newfilter.Vega_zero_flux.value*1e9, 
#                     'ab.flux.lambda':newfilter.AB_zero_flux.value*1e9, 
#                     'center_wv':newfilter.leff.value,
#                     'defualt_system':"ab"}
# )

from astropy.coordinates import SkyCoord
class Airmass:
    def __init__(self,observatory='paranal'):
        
        from astropy.coordinates import SkyCoord, EarthLocation, AltAz
        from astropy.time import Time
        from astropy import units as u
        from astroplan import Observer
        import numpy as np

        # 定义VLT的位置
        self.location = EarthLocation.of_site(observatory)  # VLT位于Paranal天文台

        # 定义观察者
        # vlt_observer = Observer(location=vlt_location, name="VLT", timezone="Chile/Continental")

    def calculate(self,
                  observation_time = Time("2024-05-23 00:00:00"),# 示例时间，可以根据需要调整
                  star_coord = SkyCoord(ra=10.684*u.deg, dec=41.269*u.deg),  # 示例坐标，可以根据需要调整

                  ):
        from astropy.coordinates import SkyCoord, EarthLocation, AltAz
        
        # 将时间和坐标转换为AltAz坐标系
        altaz = AltAz(obstime=observation_time, location=self.location)
        star_altaz = star_coord.transform_to(altaz)

        # 计算AIRMASS
        airmass = star_altaz.secz  # secz 是空质量因子

        # 打印结果
        print(f"At time {observation_time}, the airmass for the star at coordinates {star_coord.to_string('hmsdms')} is {airmass:.2f}")

        # 检查secz是否为无穷大（表示天体在地平线以下）
        if np.isinf(airmass):
            print("The star is below the horizon at the given time.")
        else:
            print(f"The airmass at the given time is {airmass:.2f}")
        return airmass.value
